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Abstract 


Analytic  gradients  of  electronic  eigenvalues  require  one  calculation  per  nuclear 
geometry,  compared  to  3n  calculations  for  finite  difference  methods,  where  n  is  the 
number  of  nuclei.  Analytic  non-adiabatic  derivative  coupling  terms,  which  are 
calculated  in  a  similar  fashion,  are  used  to  remove  non-diagonal  contributions  to  the 
kinetic  energy  operator,  leading  to  more  accurate  nuclear  dynamics  calculations  than 
thosethat  employ  the  Born-Oppenheimer  approximation,  i.e.,  that  assume  off-diagonal 
contributions  are  zero.  The  current  methods  and  underpinnings  for  calculating  both  of 
these  quantities  for  MRCI-SD  wavefunctions  in  COLUMBUS  are  reviewed.  Beforethis 
work,  these  methods  were  not  available  for  wavefunctions  of  a  relativistic  MRCI-SD 
Hamiltonian.  A  formalism  for  calculating  the  density  matrices,  analytic  gradients,  and 
analytic  derivative  coupling  terms  for  those  wavefunctions  is  presented.  The  results  of  a 
sample  calculation  using  a  Stuttgart  basis  for  K  He  are  presented.  Density  matrices 
predict  the  MRCI  eigenvalues  to  approximately  10  10  hartree.  Analytic  gradients  match 
finite  central-difference  gradients  to  within  one  percent.  The  non-adiabatic  coupling 
angle  calculated  by  integrating  the  radial  analytic  derivative  coupling  terms  matches  the 
same  angle  approximated  by  the  Werner  method  to  within  0.02  radians.  Non-adiabatic 
energy  surfaces  for  K  He  are  presented. 
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GRADIENTS  AND  NON-ADIABATIC  DERIVATIVE  COUPLING  TERMS  FOR  SPIN-ORBIT 


WAVEFUNCTIONS 


I.  Introduction 

Computational  chemistry  is  a  useful  theoretical  companion  of  spectroscopy  and 
other  experimental  approaches  to  understand  the  complex  interaction  between  atoms 
and  molecules.  It  seeks  to  explain  the  chemical  phenomena  observed  by  experimental 
science  and  further  to  predict  the  nature  of  future  experimental  results.  Over  the  past 
several  decades,  computers  have  grown  in  size  and  speed,  and  the  methods  of 
computational  chemistry  have  matured.  As  needs  are  identified,  the  field  of 
computational  chemistry  has  grown  to  meet  those  needs.  For  example,  the  coupling 
between  electronic  states  and  nuclear  degrees  of  freedom,  quantified  by  Derivative 
Coupling  Terms  (DCT)  was  largely  ignored  until  recently.  Now  that  field,  known  as  non- 
adiabatic  chemistry,  is  growing,  and  continues  to  explain  a  number  of  phenomena. 

Suppose  we  have  a  system  of  interest  which  involves  a  transition  from  one 
energy  state  to  another.  The  Diode-Pumped  Alkali  Laser  (DPAL)  is  one  such  system  in 
which  non-adiabatic  dynamics  plays  an  important  role.  In  a  DPAL,  the  alkali  atoms  are 


1 


pumped  to  the  excited  2 Pv  state  from  which  they  transition  to  the  2 Pv  state  via 

/2  /2 

collision  with  a  gas,  creating  a  population  inversion  in  the  2 Pv  state;  finally  they  lase  to 

the  ground  state  (see  schematic  in  Figure  1).  To  model  the  transition  of  the  system,  the 
nuclear  wavefunction  can  be  represented  as  a  two-component  vector, 


\ 


(i) 


and  the  Hamiltonian  which  will  propagate  it  can  be  represented  as 

h  =  k+v  =  -£— 
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where  a  indexes  the  nuclei  and  the  s's  are  the  potential  energy  values  (determined 
by  the  electronic  Schrodinger  equation).  We  restrict  the  discussion  that  follows  to  a 
diatomic  system  consisting  of  an  alkali  atom  and  a  noble  gas  atom.  Because  this 
Hamiltonian  is  completely  diagonal,  a  wavefunction  that  is  prepared  in  the  2 Py  state 


2.  Collisionally  de-excite 


1.  Pump 

3 .  Lase 

t 

2^ 


25 


Figure  1.  Schematic  of  DPALS 
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will  never  couple  to  the  2 Pv  state  (note  that  this  Hamiltonian  does  not  account  for 

/2 

interaction  with  an  electromagnetic  field  orvacuum  fluctuations).  This  Hamiltonian  is 
adiabatic  in  that  it  does  not  allow  a  wavefunction  to  transition  between  potential 
energy  surfaces,  and  is  a  result  of  the  Born-Oppenheimer  Approximation  which 
completely  decouples  the  electronic  and  nuclear  dynamics. 

Such  systems  cannot  be  studied  under  the  Born-Oppenheimer  Approximation, 
and  a  more  realistic  albeit  complicated  Hamiltonian  must  be  employed: 


(3) 


The  off-diagonal  values  in  the  kinetic  energy  operator  are  the  DCTs  which  are  a  result  of 
the  coupling  of  electronic  and  nuclear  dynamics.  They  are  not,  in  general,  negligible. 


but  they  do  allow  a  nuclear  wavefunction  to  transition  between  the  two  levels. 


Unfortunately,  these  terms  also  introduce  off-diagonal  derivative  operators  into 
the  kinetic  energy  operator,  making  the  propagation  equations  difficult  to  solve.  To 


remedy  this  shortcoming,  we  may  introduce  a  unitary  transformation,  which  is 


calculated  from  the  DCTs,  to  the  Hamiltonian  to  diagonalize  the  kinetic  energy  operator: 
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This  transformation  results  in  two  new  diabatic  potential  energy  curves,  sx  and  s2 ,  as 


well  as  a  coupling  curve,  sl2.  This  diabatic  form  of  the  Hamiltonian  still  allows 
transitions  between  the  two  states,  but  now  results  in  tractable  nuclear  dynamics 
equations. 

In  order  to  construct  this  form  of  the  Hamiltonian,  we  must  be  able  to  calculate 
the  DCTs,  Pa,  and  from  them  the  transformation  matrix  U.  Shepard's  method  of 
analytic  gradients  [1]  [2]  [3]  [4]  and  Lischka's  adaptation  of  it  to  DCTs  [5]  [6]  has  made 
their  calculation  possible  in  COLUMBUS  [7]  [8]  [9]  [10],  a  popular  quantum  chemistry 
code  package.  While  Yabushita  has  implemented  the  ability  to  calculate  the 
wavefunctions  and  energies  of  open-shell  systems,  including  spin-orbit  systems,  in 
COLUMBUS  [11],  the  analytic  calculation  of  the  DCTs  (and  energy  gradients)  for  such 
systems  has  not  been  available  for  Configuration  Interaction  (Cl)  calculations. 

Previously,  a  number  of  methods  were  available  to  approximate  inclusion  of 
spin-orbit  coupling  into  DCTs.  For  example,  code  is  available  to  calculate  spin-orbit  DCTs 
for  a  Mulit-Configu ration  Self-Consistent  Field  (MCSCF)  [12];  however,  MCSCF 
calculations  are  generally  smaller  than  Cl  calculations,  and  thus  produce  less -accurate 
wavefunctions.  Spin-orbit  coupling  can  also  be  added  perturbatively  to  the  non- 
relativistic  Hamiltonian  (see,  e.g.,  [13]  [14]  [15]);  however,  perturbative  approaches  are 
generally  not  applicable  to  larger  atoms.  It  is  also  possible  to  perform  a  spin-orbit 
calculation  to  create  the  Cl  wavefunctions,  but  then  use  a  non-relativistic  approach  to 
calculate  the  DCTs  and  resultant  transformation  which  are  applied  to  the  entire 


4 


wavefunction  [16]  [17]  [18];  however,  this  method  is  not  as  rigorous  as  including  the 


spin-orbit  effects  in  the  DCTs. 

Werner  introduced  a  rudimentary  yet  extremely  fast  way  to  calculate  the  non- 
adiabatic  coupling  angle  by  examining  the  Cl  coefficients  of  the  wavefunction,  the 
derivative  of  which  is  the  DCT  [12].  Because  of  the  speed  and  ease  of  calculation,  we 
will  use  this  last  method  as  a  quick  validation  of  the  method  presented  herein. 

The  recent  experimental  work  at  the  Air  Force  Institute  of  Technology  on  DPALs 
has  given  impetus  to  the  theoretical  effort  in  the  system's  nuclear  dynamics,  and  hence 
a  need  for  accurate,  open-shell  spin-orbit  DCTs.  This  dissertation  reviews  the 
established  methods  available  in  COLUMBUS  and  then  presents  the  formalism  whereby 
they  can  be  augmented  to  produce  analytic  gradients  and  DCTs  for  open-shell  spin-orbit 
wavefunctions.  This  formalism  has  been  implemented  in  an  experimental  version  of 
COLUMBUS  at  the  Department  of  Defense  Super-computing  Resource  Center  (DSRC)  at 
Wright-Patterson  Air  Force  Base,  OH  [19].  The  alkali-noble  gas  system  of  K  He  has  been 
chosen  to  validate  the  code  due  to  its  simplicity,  the  availability  of  a  spin-orbit  potential, 
and  its  potential  application  to  DPALs.  The  resultant  gradients,  DCTs,  and  adiabatic 
potential  energy  surfaces  of  K  Fie,  the  last  of  which  may  be  used  for  the  nuclear 
dynamics  calculation,  are  presented  herein. 
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II.  Background 


Second-Quantization  of  the  Fine-Structure  Hamiltonian 

The  energy  operator,  or  Hamiltonian,  will  be  the  principal  quantum-mechanical 
tool  used  in  this  paper  to  derive  the  quantities  of  interest  to  us,  viz,  energy  eigenvalues, 
energy  gradients,  and  derivative  coupling  terms.  The  Hamiltonian  obeys  the  eigenvalue 
equation 

//|'P/)  =  £/|T'/)  (5) 


in  which  VF/  are  the  solution  wavefunctions,  and  the  E1  are  the  corresponding  energy 
eigenvalues.  A  brief  introduction  to  the  Hamiltonian,  and  the  Schrodinger  equation,  is 
discussed  in  AppendixA. 

The  Electronic  Hamiltonian 

For  the  atomic  and  molecular  systems  which  we  will  explore  in  this  paper,  we 
can  formulate  a  Hamiltonian  in  which  the  kinetic  energy  of  each  particle  and  each  pair¬ 
wise  Coulombic  potential  energy  are  represented,  all  of  which  are  linearly  superimposed 
[20]: 


a  2ma  j  2  a./3\ra  rp\  i,j  \ri  rj\  a,i  \ra  r\ 

where  Za  and  ma  are  the  charge  and  mass,  respectively,  of  the  a"1  nucleus;  rt  and  ra 
are  the  position  of  the  ith  electron  and  ath  nucleus,  respectively;  and  V2a  and  V2  are 
the  Laplacian  operators  of  the  a,h  nucleus  and  the  ith  electron,  respectively.  In  this 
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equation  and  for  the  remainder  of  the  paper  we  will  use  atomic  units  for  which  the  mass 


of  the  electron,  the  charge  of  the  electron,  and  fi  are  all  equal  in  magnitude  to  unity. 
The  terms  in  equation  (6)  are,  from  left  to  right: 

1.  The  nuclear  kinetic  energy  summed  over  the  nuclei 

2.  The  electronic  kinetic  energy  summed  over  electrons 

3.  The  nuclear  potential  energy  summed  over  nuclear  pairs 

4.  The  electronic  potential  energy  summed  over  electron  pairs 

5.  The  nuclear-electronic  potential  energy  summed  over  nucleus-electron  pairs 
This  Hamiltonian  is  often  split  into  a  nuclear  Hamiltonian 


H 


n 


(7) 


and  an  electronic  Hamiltonian 


i  2 


ri~rj\ 


-r 


(8) 


Notably  missing  from  either  Hamiltonian  is  the  nuclear  potential  term.  This  term  will  be 
the  same  for  all  electronic  eigenvalues,  so  it  will  not  affect  the  calculations  to  come; 
rather,  it  will  be  added  to  the  final  energy  as  an  offset. 

Effective  Core  Potentials 

The  above  electronic  Hamiltonian  is  known  as  an  all-electron  Hamiltonian,  since 
it  is  assumed  that  i  indexes  all  electrons  belonging  to  the  atoms  in  the  system.  For 
larger  systems  or  small  systems  with  larger  atoms,  the  two-electron  interaction  term 
can  become  quite  cumbersome  to  calculate.  It  is  expedient,  then,  to  replace  some  of 


7 


the  core  electrons  that  play  a  lesser  role  in  the  chemistry  with  an  Effective  Core 
Potential  (ECP),  while  dealing  with  the  more  important  valence  electrons  individually 
[21]  [22],  The  resultant  electronic  Hamiltonian  is 


i  L  ‘J 


r.~rj 


■Zi 


YuECF 


00 


(9) 


where  i  indexes  the  valence  electrons  and  possibly  some  subset  of  the  core  electrons. 
The  ECP  captures  the  approximate  effects  of  the  core  electrons  on  the  valence  set. 

Relativistic  Effects 

The  above  electronic  Hamiltonian,  while  powerful,  does  not  take  into  account 
the  effects  of  relativity  on  the  energy  eigenvalues.  Appendix  B  addresses  the  Breit-Pauli 
Hamiltonian,  a  more  accurate  and  relativistic  energy  operator.  While  we  will  not 
explicitly  usethe  Breit-Pauli  Hamiltonian  in  our  calculations,  we  will  be  using  one  of  the 
quantities  derived  from  it.  We  include  the  largest  relativistic  contribution  in  the  valence 
region,  the  spin-orbit  contribution,  by  means  of  a  Relativistic  ECP  (RECP), 


H 


l  ^  ij 


ri~rj  I 


+Z uRECP (ij) 


(10) 


which  is  calculated  by  transforming  the  solutions  to  the  Dirac-Hartree-Fock  equations 
[23]  [24],  The  RECP  is  calculated  in  two  pieces:  an  Averaged  Relativistic  Electron 
Potential  (AREP),  which  approximates  the  average  fine-structure  contribution,  and  a 
Spin-Orbit  Potential  (SOP),  which  effects  the  actual  splitting  of  the  energy 
eigenfunctions.  The  spin-orbit  potential  may  be  defined  as  [11] 
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(11) 


HXo  =  'L^i(i’a)drh 

a,i,l 

in  which  cr  is  the  spin  vector  of  the  ith  electron,  lia  is  the  orbital  angular  momentum 
vector  of  the  i,b  electron  with  respect  to  the  a"' nucleus,  and  i,a )  is  a  numerically- 


calculated  radial  potential  term.  (Compare  this  term  to  the  term  hsg  in  equation  (306) 
of  Appendix  B).  Making  these  corrections  to  the  electronic  Hamiltonian  in  equation  (9) 
yields  a  spin-orbit  electronic  Hamiltonian: 


V2  1 

*.=-Z^-+*Inr 


-  + 


U  I  ri~rj\ 


\r  —  r\ 

a,i  y a  i\ 


YJu'REr(n)+H, 


(12) 


One-  and  Two-Electron  Operators 

The  electronic  Hamiltonian  in  equation  (12)  has  two  types  of  integrals  with  which 
we  will  be  concerned:  the  one-electron  operators  (those  that  sum  only  over  one 
electronic  coordinate)  and  the  two-electron  operators  (those  that  sum  over  two 
electronic  coordinates).  We  define  the  operators 


f  z  A 


\\ra  -Tiy 


+  UAREP(ri) 


hSO{i)^'Z(%i(i'a)CJrTia) 

a, l 


§(ij) 


\ri~rj\ 


(13) 


so  that  the  electronic  Hamiltonian  may  be  written  in  the  more  compact  form 


H  =  Y,h{i)  +  Yj^° 


IJ 


(14) 
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Second  Quantization 


The  Hamiltonian  introduced  in  the  previous  section  is  the  first-quantized 
Hamiltonian.  Some  of  the  distinguishing  features  of  the  first-quantization  are  [3]: 

1.  It  is  constructed  from  the  fundamental  characteristics  of  the  particles  (mass, 
charge,  position). 

2.  It  depends  on  the  number  of  electrons,  N,  and  thus  it  confines  its  eigenfunctions 
to  have  N  electrons. 

3.  There  is  no  requirement  that  the  functions  on  which  it  operates  be 
antisymmetric  or  even  built  from  an  orbital  basis. 

4.  The  choice  of  orbital  basis  has  no  effect  on  the  form  of  the  Hamiltonian. 

In  contrast,  we  wish  to  form  a  second-quantized  Hamiltonian  that  is  constructed  from  an 
orbital  basis  and  is  unconfined  by  electron  number.  Furthermore,  this  Hamiltonian  will 
be  represented  in  a  space  of  antisymmetric  many-electron  functions,  known  as 
Configuration  State  Functions  (CSFs). 

Consider  the  spin-orbital  projectors, 

lr(/)//(ryi)><r(/)//(ry;.)l  (15) 

where  r(i)  signifies  the  ith  electron  in  the  rh  spatial  orbital  and  //(«,.)  signifies  the 

same  electron  in  the  ju,b  spin  state.  Through  completeness, 

YJ\r(i)/j(coi))(r(i)/u(coi)\  =  l  (16) 

r,/u 

which  we  can  now  insert  judiciously  into  the  first -quantized  Hamiltonian: 
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^=ZZlr(0^(®i)><r(0^(®/)l,i(0Zlj(0i'(®i))<s(0v(®i)1 

i  r,/u  s,v 

+  2EElr(O^(^)><r(O^(^)l^0‘’-/)Sls(Ov(<yI)><5(Ov(^)1  (1?) 

i,j^i  r,fi  5,v 

+ Z  Z1  r  (0  /*  (®/  ))<r  (0^  (®/ )  w°  (0Z1 5  (0v  (®/  ))<s  (O' v  (®« ) 1 

i  r,/u  ^,1/ 

A  word  of  caution:  equation  (16)  is  only  true  for  a  complete  (i.e.  infinite)  set  of  orbitals. 
All  computational  chemistry  calculations  use  a  finite  basis,  and  thus  equation  (16)  is  not 
necessarily  true.  The  result  of  this  approximation  is  that  the  second-quantized 
Hamiltonian  at  the  end  of  this  section  is  not  exactly  equal  to  the  first -quantized 
Hamiltonian  at  the  beginning.  However,  in  orderto  use  a  first -quantized  Hamiltonian  in 
calculations,  it  too  must  at  some  point  be  projected  into  a  working  CSF  basis;  if  this  basis 
is  the  same  as  the  one  found  in  equation  (16),  then  the  first-  and  second-quantized 
Hamiltonians  will  havethe  same  representation  [3], 

Since  the  sums  over  spin-orbitals  are  independent,  they  can  be  combined. 
Further,  since  h(i)  and  g(i,  j)  are  not  functions  of  spin,  the  spins  are  not  involved  in 
those  integrals: 

H = Z  Z  1  r )><r(0  h (0 1 5 (0></*(®/ )lv(^ )><s(i'M®» ) 1 

i  r,n,s,v 

+  iZ  Z  \r(i)Ju((oi))(r(i)\g(<i,j)\s(i))(/u(G)i)\v(coj))(s(i)v((oi)\  (18) 

ij^i  r,/u,s,v 

+ Z  Z  1  r (0 v (®/ )><r (0 v (®/ ) hSO  (0 1 5 (O' v (®i )><5 (O' v (®» ) 1 

i  r,ju,s,v 

For  compactness,  let  us  rename  the  integrals: 

K  =  (r(i)\h{i)\s(i)) 

Srs  (j)  =  (r(0  1  8  (tj)  1  5  (0>  (19) 
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Note  that  since  g  started  as  a  two-electron  operator,  it  still  has  an  electron  dependence 


after  a  single  integration.  The  orthonormality  of  the  spin  functions  eliminates  one  of 
the  spin  sums  in  the  non-spin  terms: 

^  =  ZZK0M^))(5(0M^)K 

i  r,ju,s 

+1Z  ZIK0M®»))(s(0/*MM-j)  <20) 

i,j*i  r,ju,s 

+Z  Z  k(0^(®i))(s(0v(®i)fc 

i  r,jU,s,v 


In  order  to  remove  all  electron  dependence  from  g  ,  we  again  apply  completeness: 

H  =  X  X 1  1  Ks  + 

i  r,p,s 

i,j*i  r,ju,s  t,p,u,<j 

+Z  Z  |r(0^(^)><j(0v(®i)1^ 


(21) 

Following  the  same  prescription  as  equations  (18)  and  (20),  the  spin  functions  p  and  a 
are  removed  from  the  integral  and  one  of  the  spin  sums  is  eliminated.  We  now  define 
the  integral 


grstu  =(t(j)\grs(j)\u(j)} 


(22) 


and  we  have 

H = Z  Z 1  r  O'M®/  )><s  (0^(®/ ) 1  /7« 

i  r  ,p,s 

+iZ  Z  T,w(^^M^s(i)^M[tU)p((DMuU)p((Djy  grstu  (23) 

i,j^i  r,ju,s  t,p,u 

+Z  Z  |r(0//(^))<5(0v/(^)l/Cr 

i  r,p,s,v 


Let  us  define  the  creation  and  annihilation  operators 
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1 


a  = 

rp 


y/N  +  li 


Zlr(0^(®«)> 


r/i 


(24) 


(see  Appendix  C).  The  product  of  a  raising  and  a  lowering  operator  with  identical  indices 
will  only  have  non-zero  results  on  the  diagonal  of  an  operator  matrix. 


f 


y 

) 

A 


•Jn + 


=  ^lr(f)//(ryi))<r(f)//(ryi)l 


(25) 


We  can  now  place  these  operators  into  the  Hamiltonian: 

H  =  y  o'  a  h 

/  j  rp  sju  rs 
r,p,s 

i,j^i  r,ju,s  t,p,u 

+  y  af  d  hso 

/  j  rp  s v  rjusv 
r,ju,s,v 


Note  that  the  pair  of  operators  for  the  non-relativistic  one-electron  integral  is  spin 
preserving,  but  the  pair  for  the  spin-orbit  integral  is  not.  Note  also  that  the  g  term 
appears  to  have  four  creation/annihilation  operators;  however,  since  j  =  i  is  not 
allowed,  one  pair  of  these  operators  is  short  of  a  complete  sum.  We  can  correct  this 


omission  by  adding  and  subtracting  the  j  =  i  term: 


H  =  y  a'  a  h 

/  j  rp  sp  rs 
r,p,s 


+  2 


r,p,st,p,u  i,j 


-  KO^  (®/ ))  {s  i})P  (®i )  k  (l)p  (®i ))  (U  (l)p  (®z  ) 

+  y  a  a  hso 

/  j  rp  sv  rpsv 


(27) 
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Now  the  second  term  includes  a  complete  integral  over  the  electronic  variable  i : 


{s{i)fi{cQi)\t{i)p{(oi))  =  SstSMp  (28) 

and  the  Hamiltonian  is  simplified  as 
H  =/  a 1  a  h 

rp  sp  rs 

r,p,s 

+  4  T ^  (a]  a  a]  a  —  d^d  8 .8  )  g,  (29) 

2  \  rp  sp  tp  up  rp  up  st  pp  J  o  rstu  '  ' 

r,p,s  t,p,u 

+  'V  d]  d  hso 

/  j  rp  sv  rpsv 
r,p,s,v 


Let  us  further  define  [11]  [3] 


d'  d 

rp  sp 


p  =  f1  f  _  J7 

C  rstu  ^rs^tu  ^ ruU  st 

Z~I  /V  +  /V 

t  =a  a 

rpsv  rp  sv 


(30) 


such  that  the  final  form  of  the  second-quantized  spin-orbit  Hamiltonian  takes  the  form 


H  =  V  F  h  +1  V  p  g  +  V  S’  us° 

11  Z^ rsnrs  —  2  Z^  rstu  O  rstu  '  Z^  ^ rpsv'1  rpsv 


(31) 


r,p,s,v 


In  a  later  section  we  will  see  that  the  Ers  are  generators  of  the  Unitary  Group.  For  this 
reason,  the  second-quantized  form  of  the  Hamiltonian  will  be  central  to  the 
computational  chemistry  techniques  we  discuss  in  later  sections. 

Density  Matrices 

One  tool  available  from  second  quantization  of  the  Hamiltonian  will  be  vitally 
important  to  the  construction  of  analytic  gradients  and  DCTs  is  the  density  matrix.  Let 
|VF/)  be  an  eigenfunction  of  the  Hamiltonian  in  equation  (31).  Then  its  energy 
eigenvalue  is 
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=  Z(>i'-|£Al'i',)+i  Z  <^,|0-l»,>+  Z  ('P,lOCl'f'; 


Since  the  integral  matrices  are  not  functions  of  electronic  coordinates,  they  commute 
with  the  eigenfunctions  and  we  have 

£'=5X('f',|£i-,l'f',M  Z  Z  'CC*',  14-1'*'.)  (33) 


We  define  the  quantities 


D„=('r,|E„|'F,; 

dr»u=('V,\e\'V, 


to  be  elements  of  the  one-electron  and  two-electron  density  matrices,  respectively.  The 
spin-orbit  density  matrix,  which  would  appear  in  the  last  term  of  equation  (33),  will  be 
addressed  in  Chapter  III.  For  many  transition  properties  between  wavefunctions, 
including  DCTs,  it  is  necessary  to  use  transition  density  matrices  which  we  similarly 
define  as 


These  matrices  will  appear  often  in  the  derivations  in  Chaper  III  as  well  as  the 


appendices. 
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Gauge  Theory 


The  many-body  problem 

Systems  of  more  than  two  bodies  are,  in  general,  not  analytically  tractable  [25], 
The  sole  example  of  a  two-body  problem  is  the  hydrogen  atom  (or  exotic  variants 
thereof)  consisting  of  one  proton  and  one  electron;  even  the  basic  hydrogen  molecule 
has  no  analytic  solution.  The  energies  and  eigenfunctions  of  the  hydrogen  atom  were 
solved  analytically  long  ago;  indeed,  this  was  one  of  the  primary  reasons  for  the 
introduction  of  quantum  mechanics  [26],  The  desire  to  understand  molecules  from 
diatoms  to  proteins  and  beyond  requires  various  types  of  approximations.  This  is  the 
many-body  problem  that  is  fundamental  to  this  field,  which  forces  us  to  look  for 
computational  solutions  to  molecular  problems. 

Another  problem  that  faces  computational  chemistry  is  the  disparate  nature  of 
the  particles  involved.  Electrons  are  light  and  fast  compared  to  the  massive,  lumbering 
nuclei  (which,  for  most  calculations,  are  modeled  as  a  single  particle).  The  vastly 
different  sizes  of  these  two  particles  leads  to  two  separate  timescales  and  separates  the 
problem  of  molecular  dynamics  into  an  electronic  and  a  nuclear  piece  (see  the 
separation  of  the  Hamiltonian,  equations  (7)  and  (8)). 

The  Born-Oppenheimer  Approximation 

The  Born-Oppenheimer  approximation  attempts  to  effect  this  separation.  In  this 
section,  let  r  represent  the  collection  of  electronic  coordinates,  and  let  R  represent 
the  collection  of  nuclear  coordinates.  Assume  we  have  a  complete  set  of  orthonormal 
electronic  wavefunctions  such  that 
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(¥i  I ¥j)  =  (36) 

The  solution  to  the  Schrodinger  equation  (equation  (5))  can  be  expanded  in  this  basis: 


l'F,)  =  2Xlr,> 

i 


(37) 


which  is  called  the  Born-Oppenheimer  expansion  [27],  As  we  will  see  shortly,  these 
coefficients  are  actually  the  nuclear  wavefunctions.  Alternatively,  we  can  write 

("  \ 


‘“71 


M 

1  73 


(38) 


7  J 

in  the  matrix  representation.  The  Born-Oppenheimer  approximation  takes  advantage  of 
the  different  timescales  discussed  in  the  last  section  by  assuming  the  Hamiltonian  can 
be  partitioned  into  fast  and  slow,  or  electronic  and  nuclear,  operators.  The  Schrodinger 
equation  can  be  cast  as 


Z(»,  +  ».-£')h»i'O  =  0 

i 


(39) 


where  Hn  is  the  part  of  the  full  Hamiltonian  that  depends  on  nuclear  operators  and  He 
depends  on  electronic  operators,  as  defined  in  equations  (7)  and  (8),  respectively. 
Considerthe  integral  equation 

(O';  I  S(A+H,-£')S»I^>=0  (40> 

i 

which  is  expanded  as 

1  I  ¥i)  +  (Vj  I  H„Eb  I  Yi)~ (Vj  I  E‘Eii  I  ¥ii)  =  0  (41) 
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In  the  first  term,  the  Hamiltonian  acts  on  the  bra  and  pulls  out  an  eigenvalue,  while  the 
E/;  does  not  participate  in  the  integration,  leaving  a  delta  function.  In  the  third  term, 
neither  E1  nor  Efi  participates  in  the  integration,  leaving  another  delta  function: 

£jEv  +H^j 1  1  Vi)  ~  £l^ij  =  0  (42) 

i 


To  interpret  the  second  term,  it  is  best  to  project  it  into  coordinate  space  by  employing 
the  identity 

^\r,R)(r,R\  =  l  (43) 

r,R 

where 

3,((R)^{r,RI3„> 

V2 

-Y-*-^>(r,R\Hn\r,R)  (44) 

«  2ma 

sj(R)^(r,R\sj\r,R) 

E1  (R)^(r,R\EI  \r,R) 


Note  that  while  the  y/i  are  strictly  functions  of  electronic  coordinates,  they  are 
parameterized  by  the  nuclear  coordinates  [27],  Since  the  total  wavefunction  must  be  a 
function  of  both  nuclear  and  electronic  coordinates,  this  leaves  the  nuclear  dependence 
in  the  S/;term.  The  form  of  the  nuclear  Hamiltonian  in  coordinate  space  comes  from 
equation  (7).  (Note  that  by  including  the  kinetic  energy  operator  but  not  the  potential 
energy  operator,  we  have  inherently  prepared  for  the  coordinate  representation  since  it 
is  the  momentum  operator  that  will  be  represented  by  a  derivative  operator.  Were  we 
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to  prepare  for  the  momentum  representation,  we  might  choose  to  keep  the  potential 


energy  operator  and  not  the  kinetic.)  In  the  coordinate  representation,  equation  (42)  is 
(if)  -  S  J  -R) [S„  ( ^)^,  ( « )]  -  )  =  0  (45) 

i,a  a 

The  derivative  operator  must  be  distributed  across  the  product  in  brackets  in  the  usual 
fashion: 

VS[^(<-;*)S„(/!)]=V„-[^(r;/i)V„3„(/i)+3a(R)V^(r;/i)] 

(46) 

=  2V^,(r;*)-VA(R)+|ir,(r;R)Vj3„(R)+a„(/i)V>,(r;/i) 
which  leads  to  the  equation 

i,a  zma  (47 

+r;(r;R)r,('';«)V*SJ1(«)+r;(r;«)Vj[r,(r;«)].S„(«))-£'(«)3,,(«)  =  0 

Each  term  in  this  Schrodinger  equation  operates  on  S n(R)  or  S (/?).  Therefore,  the 
terms  can  be  collected  as  follows: 


I^v«+£, W-E'(*) 


a  1  2ma 


=,(«)= 


YJ^\dr[2¥](r,R)V^Xr,R)V^y,](r,R)Vl¥l(r’R)y-,AR) 


(48) 


J  2  ma 


The  left-hand  side  of  this  equation  is  devoid  of  electronic  dependence,  while  the  right- 
hand  side  involves  operators  which  take  derivatives  of  electronic  functions  with  respect 
to  nuclear  coordinates.  We  can  compactly  rewritethe  right-hand  side  as 

XAj,(r,R)E,(R)  (49) 

i 


where 
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(50) 


AAR)s'L4r2pAR)^.+QAR) 

a  Ama 

P“  (R)  =  j dry)(r;R)Vay,  (r.R)  o  (yt  | V„r,) 
fi"  (A)  =  |  dry)  (r.R)  V>,  (r;R)  o  |  V>,) 

leading  to  the  more  compact  Schrodinger  equation 

+  *,(*)-£'  (R)k  (R)  =  £>„  (R)3„  (*)  (51) 

L  a  2ma  J 

or,  in  matrix  representation, 

-X-i-IV=+e(R)-I£'(fi)  S,(R)  =  A(R)S,(R)  (52) 

«  2ma 

Note  that  all  matrices  on  the  left  are  diagonal  matrices,  while  the  matrix  on  the  right 
side  is  unrestricted.  If  we  assume  that  derivatives  of  electronic  functions  with  respect  to 
nuclear  coordinates,  or  DCTs,  are  negligible  compared  to  the  rest  of  the  equation,  the 
right  side  of  equation  (51)  or  (52)  is  summarily  zero,  and  we  have  the  adiabatic  nuclear 
wave  equation  [28]: 

-E^-IV»+£W-K'W  S,(R)  =  0  (53) 

«  2  ma 

which  is  the  extent  of  the  Born-Oppenheimer  approximation;  nuclear  and  electronic 
dynamics  have  been  completely  separated. 

Adiabatic  chemistry,  the  subject  of  most  computational  chemistry,  relies  on  the 
Born-Oppenheimer  approximation  and  the  resultant  adiabatic  equation  above.  The 
lumbering  nuclei  are  frozen  in  space  and  the  electronic  Schrodinger  equation  is  solved 
around  them.  The  eigenvalues  of  that  equation,  £,  parameterized  by  the  nuclear 
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coordinates,  serve  as  the  adiabatic  potential  energy  surface  in  equation  (53).  (The  term 


surface  is  used  loosely,  since  it  will  have  as  many  dimensions  as  there  are  internal 
nuclear  coordinates.)  This  approximation  is  feasible  because  the  timescale  of  the 
electrons  is  so  short  that  the  nuclei  are  shepherded  by  the  electronic  wavefunction 
(rather  than  by  information  about  the  precise  location  of  individual  particles),  which 
adjusts,  for  practical  calculations,  instantaneously  to  the  nuclear  motion.  Since  all 
operators  in  the  nuclear  equation  are  diagonal,  these  equations  are  simply  solved,  and 
the  nuclear  dynamics  can  be  studied  with  relative  ease. 

Non-adiabatic  Chemistry 

Though  adiabatic  chemistry  provides  a  straightforward  way  to  tackle  nuclear 
dynamics,  the  Born-Oppenheimer  approximation  is  a  compromise.  The  term  adiabatic 
in  this  sense  refers  to  the  fact  that  the  complete  separation  of  fast  and  slow  variables 
does  not  allow  the  nuclear  wavefunction  to  pass  from  one  potential  energy  surface  (or 
set  of  electronic  eigenvalues)  to  another  quantum  mechanically,  a  diabatic  process,  by 
transferring  electronic  potential  energy  to  nuclear  kinetic  energy.  This  means  that  a 
molecule  has  no  way,  apart  from  electromagnetic  interaction,  of  changing  its  electronic 
state.  We  know  this  premise  to  be  false,  as  chemistry  is  wrought  with  examples  of 
molecules  which,  after  briefly  contacting  another,  rebound  with  a  higher  (or  lower) 
energy  level.  For  problems  in  which  these  transitions  are  critical,  the  field  of  non  - 
adiabatic  chemistry  is  more  appropriate. 

Non-adiabatic  chemistry  still  relies,  at  least  partially,  on  the  Born-Oppenheimer 
approximation.  It  still  requires  that  the  Hamiltonian  be  separated  into  nuclear  and 
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electronic  operators,  and  that  the  wavefunction  be  a  product  of  nuclear  and  electronic 


wavefunctions  (this  latter  requirement  being  fundamental  to  gauge  theory,  as  explained 
shortly);  however,  it  relaxes  the  assumption  that  the  elements  of  the  DCT  matrix  in 
equation  (52)  are  negligible.  To  understand  where  these  terms  are  not  negligible, 
considerthe  following  derivation,  beginning  with  the  electronic  Schrodinger  equation 
[29]: 

He  I  if/,)  =  e,  \y/l).  (54) 

Taking  the  nuclear  derivative  of  this  equation  and  integrating  against  y/]  where  j^i 
yields 

(Yj  I  V//  I  y/t)  +  (i//j\HeW y/t)  =  { y/.  Wei  I  y/.)  +  {y/j  I  st  I V y.)  (55) 

where 

V-2X  (56) 

a 

In  the  second  term,  the  electronic  Hamiltonian  acts  on  the  bra  pulling  out  an  eigenvalue; 
in  the  third  term  the  derivative  of  the  eigenvalue  is  not  involved  in  the  integration,  and 
the  integral  evaluates  to  zero  by  orthogonality;  in  the  fourth  term  the  eigenvalue  is  not 
involved  in  the  integral,  but  it  does  not  evaluate  to  zero: 

(Yj  I V//,  I  y/,)  +£j(Yj  I V  Yi)  =  £,(Yj  I V (//,)  (57) 

The  result  gives  an  alternate  definition  to  the  elements  of  the  P  matrix  part  of  the  DCT, 
which  were  defined  in  equation  (50): 
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(58) 


where 


P,  =  (y/j  I V  i//t)  = 


(Yi\  VHe\V/i) 


8,  - 


a 


a 


(59) 


(Note  that  P  is  a  vector;  with  this  understanding,  we  will  drop  the  arrow  notation). 
Thus  we  see  that  elements  of  P  and  hence  A  become  significant  when  the  energy 
eigenvalues  of  two  distinct  wavefunctions  become  very  close.  Note  that  when  si  =  s  jt 
P..  is  undefined;  such  points  are  called  conical  intersections  [30],  For  the  diatomic 
systems  that  will  be  the  focus  of  this  paper,  such  crossings  occur  only  between 
wavefunctions  of  dissimilar  spatial  symmetries,  which  will  not  be  subject  to  derivative 
coupling)  [31]. 

Derivative  Coupling  Terms 

Given  that  there  are  certain  points  in  the  nuclear  manifold  where  the  DCTs 
become  significant,  they  deserve  closer  scrutiny.  Consider  again  the  matrix  equation 
(52),  now  including  the  definitions  (50): 


IVl+z-lE 


IV« ) + Q“ 


(60) 


where  the  nuclear  coordinate  dependence  of  operators  and  functions  is  understood  and 
has  been  dropped  from  the  notation  for  simplicity.  Since  the  DCT  has  second  nuclear 
derivatives  in  it,  it  is  part  of  the  kinetic  energy  term  rather  than  the  potential: 
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In  this  form  of  the  equation  both  P“  and  Q“  have  non-zero  off-diagonal  derivative 
components  which  make  this  set  of  equations  very  difficult  to  solve  in  its  current  form. 
We  require  a  similarity  transformation  which  will  diagonalize  the  kinetic  energy  term 
but  may  undiagonalize  the  potential  energy  term  (g) ,  a  compromise  which  yields  an 
equation  with  no  off-diagonal  derivative  operators.  This  transformation  takes  the 
equation  from  the  adiabatic  basis  to  a  diabatic  basis.  Diagonal  elements  of  £  in  this 
diabatic  basis  are  the  diabatic  potential  energy  surfaces  and  are  coupled  by  the  off 
diagonal  elements  of  £  .  This  transformation  from  the  adiabatic  to  the  diabatic  basis  is 
called  the  non-adiabatic  transformation. 

In  order  to  find  that  transformation,  it  is  first  necessary  to  manipulate  the  form 
of  the  DCTs  [27],  Let  us  first  derive  an  alternate  form  of  Q“ .  In  order  to  accomplish 
this,  we  first  calculate  the  divergence  of  P“  : 

(IV«-P *\=Va^dry,]{r)Vay,i{r) 

=  J drWa  (y/)  (r))-  Vay/i  (r)  +  ^  dry/]  (r)  V2ay/i  (r)  (62) 

=  \drVa(y/*(r))-vay/i(r)  +  Q“ 

Now  consider  P“2 : 

(p“  ‘ p" )..  =  Xj dr'¥  ir')Wa¥k  (''Oj drwl  {r)^a¥j  (r)  (63) 

1  k 

Since  P“is  antisymmetric, 

K  =-K  ->  jdr^y,^  =-\drVya<ri  (64) 
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and  so  it  follows  that 


(pa  '  P“  ),  =  ~Y,\dr  'W a  \_Vi  (r')-  J dry/’  (. r)Way/j  (r)  (65) 

1  k 

Furthermore,  if  the  electronic  eigenfunctions  are  real, 

(p“  •?“)..  =  ~Yj\dr'V »  \_¥i  (r')\Vk  (r  ')J drwl  {r)Way/j  (r)  (66) 

1  k 

which  can  be  rewritten  as 

(p«.p«)  f  jrfjr,^[^(r,)^(r)]va[^;(r')]va^.(r)  (67) 

1  k 

The  double  integral  over  the  sum  of  k  -indexed  functions  is  a  completeness  identity 
(analogous  to  ¥k)(Vk  I  )>  and  produces  a  dirac  delta  function: 

k 

-J dr'Y\yk  (/■  V*  (/■)]  Va  \y/*  (/■')]' SI ay/j  (r)  =  - \dr\dr' S(r-r') Va  \y.  (r')] Vay/j  (r)  (68) 

k 

and  thus  we  conclude  that 

(P"'P")„  =-J<frV„r;(TV„r,(r)  I69) 

Thus -(P" -P“)  is  the  term  added  to  Q"  in  equation  (62).  This  means  that  Q“  and 

V  )ij  J 

P“  are  related  by 

Q“  =  (lVa  •  P“  )  +  P“2  (70) 

We  can  usethis  relationship  to  eliminate  Q“  in  equation  (61): 

-V— —  TlV2  +2(P“-IV  )  +  (lV  •Pa')  +  Por21  +  s  — I E1  S,=0  (71) 

a  2maL  J 
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Since  Pa  and  Vffare  non-commuting  operators,(P“ -P") .  Furthermore, 
note  that 

IV„  ■(P"H,)  =  (lV„  -P")S,  +P"  -IV„E,  (72) 

so  that  the  operator  relationship 

(lVg  •  P“ )  =  IV a  P“  -  P“  TVa  (73) 

can  be  substituted  into  equation  (71),  yielding 


Note  that  this  equation  looks  like  the  adiabatic  nuclear  Schrodinger  equation 
(equation(53)),  with  the  exception  that  -f^IVa  +P“  has  replaced  -i  IP.  as  the 

a  a 

momentum  operator.  This  form  is  called  the  gauge-covariant  momentum  [32]  for 
reasons  that  will  be  made  clear  in  the  next  section  (cf.  the  momentum  operator  in  the 
Dirac  equation,  equation  (289)  in  Appendix  B). 

This  simplification  has  not  changed  the  difficult  nature  of  the  equations  as  P“ 
still  has  off-diagonal  derivative  operators.  A  transformation  is  still  necessary;  however, 
with  the  Pa  identified  in  this  form  as  the  gauge  potential  [32](in  electrodynamics  the 
vector  potential  plays  the  same  role),  we  can  use  the  constructs  of  gauge  theory  to  find 
that  transformation. 
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Derivative  Coupling  Terms  in  Gauge  Theory 


By  separating  the  wavefunction  and  Hamiltonian  into  nuclear  and  electronic 
parts,  the  Born-Oppenheimer  approximation  has  set  the  stage  for  casting  non-adiabatic 
chemistry  in  the  language  of  gauge  theory.  The  nuclear  coordinates  serve  as  an  external 
manifold,  which  is  a  Hilbert  space.  At  each  point  on  the  manifold  there  is  an  internal 
space,  which  is  the  space  spanned  by  the  electronic  wavefunctions.  Hence  the  point  on 
the  nuclear  manifold  parameterizes  the  electronic  wavefunctions.  Together,  these 
spaces  form  a  fiber  bundle,  for  which  the  derivative  coupling  matrix  serves  as  the  gauge 
potential.  If  the  internal-space  basis  did  not  change  as  a  function  of  nuclear 
coordinates,  the  gauge  potential  would  be  everywhere  zero,  and  the  fiber  bundle  would 
be  considered  trivial.  This  is  the  algebraic  equivalent  of  the  complete  Born- 
Oppenheimer  approximation  which  considers  those  changes  negligible.  But,  as  the 
electronic  wavefunctions  do  change  with  nuclear  position,  the  fiber  bundle  is  non-trivial, 
and  the  gauge  potential  serves  as  a  connection  that  dictates  how  the  internal  space 
slowly  rotates  as  a  function  of  external  coordinates.  According  to  fiber  bundle  theory, 
the  bundle  at  each  point  on  the  manifold  is  equivalent  to  any  other,  and  they  can  be 
transformed  one  into  the  other  by  a  gauge  transformation.  That  is,  although  the 
electronic  basis  functions  at  different  nuclear  positions  are  not  the  same,  they  differ 
only  by  a  unitary  transformation,  which  will  be  the  non-adiabatic  transformation 
mentioned  previously.  By  understanding  the  underlying  theory  of  gauge  potentials,  we 
can  see  the  utility  of  the  derivative  coupling  terms  in  effecting  that  transformation . 
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Gauge-covariant  derivative 


In  this  section  we  introduce  an  alternative  perspective  on  the  gauge  covariant 
derivative  in  equation  (74).  The  total  differential  of  the  vector  in  equation  (38), 


produces  a  vector  representing  the  change  in  the  components  of  H/ ;  however,  this  is 
not  the  total  change  of  |XF/),  as  it  does  not  account  for  a  change  in  the  basis  [32],  The 
total  change,  represented  by  the  total  differential  of  the  wavefunction,  using 
equation(37),  is 

d  I  ¥/>  =  Yj(dEn  1  Vi)  +En  1  dVt))  (76) 


where  the  Ez  are  the  nuclear  wave  functions  and  the  | y/^  are  the  electronic 
wavefunctions,  and  the  total  differentials  are  defined  as 

d  ''*'/> 3 

a.k  GKak 

d 

d'R,  =  /  - ’R.clR  , 

Ii  a  p  Ii  ak 

a,k  CKak 


(77) 


The  first  term  under  the  summation  in  equation  (76)  looks  like  a  standard  differential, 


i.e.,  it  takes  the  differential  of  each  component  of  the  eigenvector  as  in  equation  (75). 
The  second  term,  however,  takes  into  account  the  change  in  the  internal  (electronic) 
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space  due  to  a  change  in  external  (nuclear)  space.  To  develop  the  meaning  of  I  dy/^ , 
consider  that  we  effect  a  change  in  I  <//.)  via  a  unitary  transformation: 

U\^i)  =  YjCj  \Wj)  (78) 

j 

Since  the  group  of  unitary  transformations  of  degree  n ,  U  (n) ,  (which  can  be 
represented  as  theset  of  all  nxn  unitary  matrices)  is  a  compact  Lie  group,  it  can  be 
expressed  as  the  exponentiation  of  the  underlying  Lie  algebra  u(n )  [33]: 

u(e1(*),«!(*)....)=n«p(^(jf)*“)  ,79) 

a 

where  the^are  the  parameters  and  the  Xfl  are  the  generators.  An  infinitesimal 
rotation  of  the  wavefunction  is  effected  by  [32] 

U (dQx  (/?),  d#2  (/?),...)  I  y/j)  =  ]^[exp(d#fl  (/?)  Xa )... I  ^,)  =1  yc{)+ 1  dy/^  (80) 

a 

This  can  be  shown  by  using  the  Taylor  series  expansion  of  exp  and  keeping  only  linear 
terms: 


t(dex(R),de2(R),..)\¥i)=  i+5>0a(*)xa  W,.)=i^>+i^>  (si) 


V  Cl 


This  implies  that 


V  a 


(82) 


We  can  now  substitute  the  left-hand  side  ofthis  equation  into  equation  (76): 
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(83) 


d\x¥I)  =  2>3»'  ^>+2^(«)xaE3«|t"-> 

i  a  i 

->rfS,+2>0„(K)X„S, 

a 


where,  in  the  second  line,  we  have  collected  the  coefficients  into  vectors  and  consider 
the  matrix  form  of  the  generators.  The  second  term  has  now  become  an  operator 
simply  multiplied  by  the  original  wavefunction.  We  recognize  that  the  total  differential 
of  the  eigenfunction  I VF/)  is  more  than  just  the  total  differential  of  the  vector 
components,  and  give  it  the  symbol  d ' : 


d\x¥I)  = 2>s„  i  Wl) +2>„  («)x,,£s„  i  r,>  -» 

i  a  i 

a 


(84) 


where  again  in  the  second  line  we  have  moved  from  the  abstract  operator-ket  form  of 
the  equation  to  the  matrix-vector  form.  Since  we  can  choose  the  nuclear  coordinates  to 
be  independent  of  each  other,  the  total  derivatives  are  equal  to  their  partial  derivatives: 


8R„, 


7  _8_ 

"I  8R„t 


+Z  sH,(k)xA 


(85) 


which  we  can  collect  into  vector  form: 


+Zv<A(*)xA 


and  rename 


a 


(86) 


(87) 


to  find  the  relationship 


\  a  J 


(88) 
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This  is  the  gauge-covariant  gradient  as  found  in  equation  (74). 


This  derivation  in  gauge  theory  reveals  several  characteristics  of  the  derivative 
coupling  terms.  First  of  all,  it  recognizes  P“  as  the  gauge  potential.  Specifically  from 
equation  (87),  we  see  that  the  matrix  of  derivative  coupling  terms  is  the  sum  of  the 
gradients  of  the  rotation  parameters  multiplied  by  their  respective  generators.  Thus  the 
DCTs  are  directly  connected  to  the  rotations  between  electronic  eigenfunctions  at 
different  points  in  nuclear  space.  Since  unitary  generators  are  antihermitian,  this 
relationship  shows  that  P“  is  antihermitian  (or  antisymmetric  in  the  case  of  real 
wavefunctions)  as  well.  Its  antisymmetry  implies  that  the  diagonal  elements,  i.e. 

(i//i  I  Va¥i),  are  zero,  which  further  implies  the  gradient  of  a  real  wavefunction  has  no 
projection  onto  the  wavefunction  itself. 

The  non-adiabatic  transformation 

The  purpose  of  our  introduction  to  gauge  theory  was  to  find  the  non-adiabatic 
transformation  to  make  equation  (74)  more  easily  solvable,  i.e.,  to  transform  it  into  an 
equation  in  which  P“  is  diagonal  (since  P“  is  antisymmetric,  its  diagonal  form  will 
actually  be  the  zero  matrix).  We  can  use  the  form  of  P“  in  terms  of  rotations  (equation 
(87))  to  find  that  transformation.  The  transformation  is  position  dependent,  and  will 
vary  based  on  nuclear  coordinates.  For  that  reason  it  is  categorized  as  a  local  gauge 
transformation  [21].  The  effect  of  applying  this  transformation  everywhere  on  the 
nuclear  manifold  will  be  to  align  all  the  internal  spaces;  that  is,  it  will  remove  the  nuclear 
coordinate  parameterization  from  the  electronic  wavefunctions.  There  will  still  be  an 
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arbitrary  but  constant  rotation  common  to  all  the  functions,  similar  to  the  choice  of  a 


potential  energy  offset.  Since  this  part  of  the  transformation  is  common  to  the  entire 
manifold,  it  is  a  global  gauge  transformation  [27]  and  can  beset  arbitrarily. 

Consider  a  global  gauge  transformation  Ut on  the  vector  IVaH/  (R)  [27]: 

(R)  =  U*IV„  (UU*E,  (R))  =  IV,  (§,  (R))  (89) 

where  UrH;  (i?)  =  E ,  ( R ) .  The  nature  of  the  del  operator  has  not  changed,  and  so  it  is 
globally  gauge  invariant.  Now  consider  a  local  gauge  transformation,  U’  ( R ) : 

U*  ( R)IV„E,  (R)=Ut  (R)IV„  (U(R)U’  (R)=,  (R))  = 

(yuj 

U,W<Va(U(R))S,(R)  +  IV„(E,(R)) 

Note  that  Ut  (i?)^IVaS/  (R)j  ^  IVa  ( (i?)),  and  so  a  gauge  transformation  changes 

the  nature  of  the  del  operator,  and  hencethe  nature  ofthe  Schrodinger  equation.  Thus 
the  bare  del  operator  is  not  gauge-covariant.  Now  consider  instead  the  transformation 

ofthe  vector  +P"  (/?))S7  (i?) : 

U*(R)(iv„+P“(R))e,(R)  = 

uf  (R)IV,  (U(R))E,  (R)+IV„  (E,  (R))  +  U>  (R)P“  (R)U(R)(S,  (R)) 

If  we  define 


P"  (R)=  Uf  (i?)P"  (i?)U(i?)  +  Uf  (i?)lVff  (U(i?))  (92) 


then  we  see  that 


U*  (R)(IV,  +P"  (R))E,  (R)  =  (IV,  +P"  (R))s,  (R)  (93) 
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and  the  form  of  the  gauge-covariant  gradient  operator  is  unchanged  by  the 
transformation.  If  we  apply  the  above  local  gauge  transformation  to  equation  (74)  we 
find 

a  z-,na 

and  indeed,  the  nature  of  the  equation  has  not  changed.  Forthese  equations  to  be  put 
into  a  form  in  which  they  are  most  easily  solved,  we  require  that  P  =  0 .  From  equation 
(92),  this  condition  implies  that 

P“(/?)U(/?)  =  -IVa(U(/?))  (95) 

Combining  this  equation  with  the  form  of  P  from  (87)  leads  to  the  equation 

2v„(?;,(«)x;„.u(r)  =  -2iv„(u(/;))  (96) 

a,m'  a,m 

Since  U ( R )  is  unitary  it  can  be  expressed  as  a  sum  of  exponentiated  parameter- 
generator  products  [34], 

U(/f)  =  2>pK(fl)X„]  (97) 

m 

and  equation  (96)  can  be  rewritten  as 

E  V  A  (J?)5C£exp[0„  (J?)X„]  =  -£IV„  exp[0„,  (*)X„]  (98) 

a,m'  m  a,m 

If  P“  is  calculated,  the  are  known  parameters  used  to  solve  for  the  0m(R) .  The 

del  operator  on  the  right-hand  side  of  the  equation  pulls  out  the  derivatives  of  the 
arguments  of  the  exponential  function  yielding 
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(99) 


S  VA  ( ■ « ) XL-  S' =*p  [0.  ( ■ R) X.  ]  =  -S' v„  K  ( ■ R)]  X.  exp \em  (r)  x„  ] 

<r,ra'  m  a,m 

a,m  a,m 

While  it  is  true  that  the  set  X'm,  and  the  set  X"  span  the  same  space,  they  are  not 
necessarily  the  same  set,  and  hence  we  cannot  conclude  that  the  parameters  &m  and  9"m 
are  the  same  for  all  m .  In  this  case,  there  is  no  clear  path  to  discover  the  parameters 
which  will  transform  P“  to  0,  and  an  approximate  method  must  be  used  [27],  The 
exception  to  this  case  is  when  only  two  surfaces  are  coupled,  and  hence  only  one 
generator  and  parameter  are  present,and  we  can  indeed  conclude  that 

6"  =  -O'  (100) 

is  the  condition  for  which  the  unitary  gaugetransformation  will  yield  P“  =0  for  all  a  . 

Consider  the  simple  example  of  coupling  a  two-level  system  with  only  one 
nuclear  coordinate.  The  DCT  matrix  is 
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( R0  is  chosen  arbitrarily,  and  determines  the  global  gauge.)  In  this  manner,  the  local 
gauge  transformations  can  be  determined  on  the  nuclear  manifold,  leading  to  solvable 
equations  of  the  form 


-Y  —  IVl+E-lE1 
„  2  in  “ 


H7=0 


(104) 


where  IV^  and  IE1  are  diagonal,  but  £  is  not.  It  is  to  this  end  that  finding  the 
derivative  coupling  terms,  i.e.  (y/i  Wi/Zj) ,  is  important. 


Wavefunctions  and  the  Graphical  Unitary  Group  Approach 

In  order  to  perform  a  calculation  of  molecular  energy  properties,  we  must  find  a 
suitable  representation  of  the  Hamiltonian.  Let  us  assume  that  a  suitable  set  of  one- 
electron  molecular  orbital  basis  functions  (MOs)  has  been  furnished  by  an  established 
Multi-Configuration  Self  Consistent  Field  (MCSCF)  procedure  [20]  [35],  From  these  MOs, 
we  produce  a  basis  of  CSFs,  which  are  antisymmetric  many-electron  functions,  in  which 
the  electronic  Hamiltonian  and  its  wavefunctions  will  be  represented.  Our  basis  of 
choice,  which  is  the  one  used  in  the  COLUMBUS  programs,  is  the  Gel'fand-Tsetlin  basis 
[7]  [8]  [9]  [10]  [36],  The  procedure  to  build  this  basis  is  known  as  the  Unitary  Group 
Approach  (UGA).  This  section  will  present  the  construction  of  the  Gel'fand-Tsetlin  basis 
via  UGA,  its  value,  and  how  it  is  implemented. 
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The  Unitary  Group 


Irreducible  Representation 

Many  of  the  characteristics  of  a  group  are  discoverable  in  abstract  or  operator 
form;  however,  often  it  is  necessary  to  use  a  matrix  representation  of  the  group  to 
better  understand  its  nature.  For  any  given  group  there  are  several  matrix 
representations.  For  example,  the  group  SU(2)  has  a  defining  representation  on  2x2 
matrices,  but  has  a  regular  representation  in  3x3  matrices  [34], 

The  most  useful  matrix  representation  we  will  encounter  for  our  purposes  will  be 
the  irreducible  representation,  or  irrep.  For  a  non-abelian  group,  not  all  elements 
commute  and  hence  not  all  can  be  diagonal  matrices  in  the  same  representation.  The 
irreducible  representation  is  the  basis  in  which  all  group  operators  are  block- 
diagonalized.  While  not  as  advantageous  as  a  purely  diagonal  representation,  the  block- 
diagonal  form  does  have  the  advantage  of  reducing  the  effort  of  nxn  matrix 
diagonalizations  to  separate lxl,  2x2,  3x3  etc.  diagonalizations  (n.b.  these  separate 
blocks  are  also  called  irreps  [37]). 

Building  the  CSF  space 

For  a  given  molecular  system,  we  assume  a  collection  of  N  indistinguishable 
electrons  and,  for  each,  a  vector  space  V  of  n  orthonormal  orbitals.  The  space  in  which 
such  a  system  exists  is  the  TV'7' -rank  tensor  product  of  V,  V®N ,  which  is  nN 
dimensional  [38],  This  is  known  as  the  carrier  space,  and  is  conceptualized  in  Figure  2. 
Given  these  two  groups  of  objects,  i.e.  electrons  and  orbitals,  there  are  two  distinct 
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ways  to  enumerate  our  basis.  In  the  first  place,  we  could  leave  the  orbitals  fixed  and 


permute  the  electrons  among  them.  In  the  second  place,  we  could  leave  the  electrons 
fixed  and  rotate  the  orbitals  around  them.  The  latter  method  is  the  Unitary  Group 
Approach,  while  the  former  is  the  Symmetric  Group  Approach  (SGA)  [39],  A  discussion 
of  the  SGA,  not  as  relevant  to  the  path  of  this  discussion,  can  be  found  in  Appendix  J. 
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Figure  2.  Visualization  of  carrier  space  of  N  electrons  and  n  orbitals 


Generators  of  the  Unitary  Group 

The  Unitary  group  is  of  particular  importance  to  the  field  of  Quantum  Chemistry 
because  of  the  manner  in  which  the  second-quantized  Hamiltonian  is  constructed.  The 
anticommutation  relation  [40]  of  the  raising  and  lowering  operators  from  definition  (24) 
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(105) 


( a*  , a  >  =  a*  a  +  d  d'  =1 

(  r/u  ’  r/u )  r/u  r/u  r/u  r/u 

(see  Appendix  C)  lead  to  the  following  commutation  relation  between  the  operators 
introduced  in  equation  (30)  [34]: 


y  a*  ,a  „y 

/  j  r/u  s/u  ’  /  7 


a!  a 

t/U  U/U 


Ed1  ,  a  ,,a]  a  =  /  ( 

_  r/u  s/u  ’  t/u  u/i  J  \ 


t  ^ 


,t  ^ 


a  ,a  ,, a,  a  \  =  >  \a  ,a  ,a  a  —a,  a  a  ,a 

ru  su  ’  tu  uu  I  /  -t  l  r/u  s/u  t/u  u/t  t/u  u/l  r/u  sju 


0 


H/1 


HH 


y  la'  ,( 8 .8  ,  —a)  a  . )  a  —a)  (d  8  ,  —  a'  a  )  a  , } 

I  r/u  \  st  ju  /u  t/u  s/u  J  u/u  tju  \  ur  ju  ju  r/u  u/u )  s/u  \ 

/u'/u 

y  (a^  , a  8 ,8  ,  —o'  ,a]  a  ,a  —a]  a  ,8  8,  +a]  a'  , a  a  , ) 

rlL  st  Ml1  rfl  tl1  SM  u/*  tf*  SM  ur  /u  /u  t/u  r/u  u/u  s/u  ) 


(106) 


/in 


=  £„A-e,A,  +  Z( 


—a  ,a,  a  ,a  +a '  a1  ,a  a 

r/u  t/u  s/u  u/u  t/u  r/u  u/u  s/u 


) 


Since  raising  operators  commute  with  each  other,  and  lowering  operators  do  as  well, 
the  terms  left  under  the  summation  sum  to  zero,  leaving 


\E  ,E,  ]  =  E  8, -E,  8  (107) 

L  rs  ’  tu  J  ru  st  ts  ur  '  • 

We  will  show  that  these  operators  represent  generators  of  the  Unitary  group. 

Let  U  be  an  element  of  the  group  U(n)  (the  set  of  all  nxn  unitary  matrices). 
Unitary  matrices  obey  the  bilinear  condition 

UU  =  1  (108) 

If  we  let 


U  =  exp(K) 

Then  we  have  the  condition  from  equation  (108) 

exp(K)t  exp(K)  =  exp(Kf)exp(K)  =  1 


(109) 


(110) 


38 


Since  U  and  U1  commute,  so  do  K  and  K1,  and 

exp  exp  (K  )  =  exp(K  +  Kf)  (111) 

which  leads  to  the  condition 

UfU  =  1— »K  +  Kt=0  (112) 

Element-wise,  this  condition  indicates  that 

Ky+(K.,)*  =  0  (113) 

or  that  K  is  antihermitian.  In  the  most  general  case  K  will  be  complex,  and  can  be  split 
into  an  antisymmetric  real  part  and  a  symmetric  imaginary  part: 

K  =  A  +  /S  (114) 

2  2  . 

VI  —  VI  VI  i  VI 

Without  further  restriction,  this  indicates  -  parameters  for  A  and - 

2  2 

parameters  for  S,  or  a  total  of  n2  parameters  for  K  .  This  is  the  same  number  of 
parameters  as  GL(n,R) ,  the  group  of  all  nxn  non-singular  real  matrices  [34], 
Consequently,  the  generators  of  U  ( n )  are  linear  combinations  of  the  generators  of 
GL(n,R) ,  and  we  can  use  either  set  of  generators  to  construct  U  .  In  the  most  general 
formulation,  a  unitary  operator  can  be  constructed  from  generators  as 

U  =  expfx^E,,l  (115) 

V  re  J 

(Cf.  equation  (79);  these  formulations  are  not  equivalent;  that  is,  the  same  parameters 
used  in  either  case  will  lead  to  different  results.  However,  they  are  both  equally  valid .) 
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Since  the  argument  of  the  exponent  must  be  antihermitian,  one  of  the  following 
restrictions  must  apply: 

1.  The  Era  are  antihermitian  (i.e.,  generators  of  U  (n) ),  while  the#ri  are 
unrestricted  real 

2.  The  Era  are  split  into  a  set  of  symmetric  and  a  set  of  antisymmetric  matrices;  the 
former  take  imaginary  parameters  and  the  latter  take  real  parameters 

3.  The  Era  are  unrestricted  (i.e.  generators  of  GL(n,R ) )  while  0rs  =  —0*r 

We  shall  choose  to  use  the  GL(n,R )  generators  which  are  unrestricted;  thus  we  will 
confine  ourselves  to  condition  3.  The  generators  Etj  of  GL(n,R )  are  sparse  matrices 
with  only  one  non-zero  entry;  specifically, 

(£■,)„  =4<5,  (116) 

(Consequently,  the  generators  associated  with  A  would  have  the  form  and 

V 

those  associated  with  S  would  have  the  form  E..+E...)  For  example,  in  matrix  form, 
the  generator  E^3  would  be 


^0  0  0 

0  0  1  ••• 

El3~  0  0  0 

v :  :  :  '■) 

It  is  clear  through  matrix  multiplication  in  this  representation  that 


EE  -  E  S 

^ij^kl  L-iluik 


and  hence  the  generators  obey  the  commutation  relation 

[E0^Eu~\  =  EaSjk-EjkSa 


(117) 


(118) 


(119) 
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This  is  the  exact  commutation  relation  we  saw  in  equation  (106),  and  thus  we  see  that 


the  second-quantized  Hamiltonian  applied  to  a  system  with  n  orbitals  is  constructed 
from  the  generators  of  U  (n) .  This  connection  is  why  the  unitary  group  will  be  so 

important  in  outlining  a  basis  in  which  to  express  the  Hamiltonian;  that  process  is  the 
UGA  [38],  Note  that  the  generators  in  equation  (30)  were  summed  over  spin,  and  so  the 
unitary  generators  do  not  distinguish  spin-up  from  spin-down.  We  can  also  create  the 
operators  [38] 

EMV  =  (120) 

r 

for  functions  of  spin  ratherthan  for  functions  of  electronic  coordinate.  They  obey  the 
same  commutation  relation;  however,  since  /u  and  v  can  be  selected  from  two  spin 
states,  these  operators  correspond  to  the  generators  of  U  (2) .  The  complete  space  in 
which  we  represent  the  Hamiltonian  for  n  orbitals  is  the  direct  product  space  of  the 
spatial  functions  and  the  spin  functions,  U (n)®U (2),  which  is  a  subgroup  of  U (2 n) . 

Note  that  the  spin  index  necessarily  remains  in  the  generators  from  which  the 
spin-orbit  Hamiltonian  is  constructed  (see  equations  (29),  (30)).  The  spin-orbit 

Hamiltonian  is  constructed  out  of  generators  from  U  (2 n),  and  the  CSFs  span  an  irrep  of 
a  subgroup  of  this  space. 

Young  frames 

Before  discussing  the  representation  of  the  unitary  group,  it  is  n  ecessary  to  make 
a  short  aside  to  the  subject  of  Young  frames,  which  will  provide  a  useful  labeling 
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scheme.  Every  integer  N  can  be  partitioned  a  number  of  ways.  For  example,  3  can  be 


partitioned  into  3,  2+1,  or  1+1+1.  Table  1  enumerates  the  partitions  of  the  first  five 
integers.  We  can  also  use  boxes  to  indicate  the  partitions,  as  shown  in  Figure  3.  These 
box  figures  are  called  Young  Frames  [39],  To  prevent  redundancy  when  creating  these 
frames,  it  is  necessary  to  impose  the  rule  that  a  row  be  equal  to  or  shorter  than  the  row 
above  it,  and  each  row  be  left -justified.  Figure  4  is  an  example  of  a  Young  frame  which 
is  not  allowed. 

These  frames  were  first  used  to  label  irreps  of  the  symmetric  group;  however, 
because  of  their  close  connection,  they  can  also  be  used  to  label  the  irreps  of  the 
unitary  group,  as  the  succeeding  discussion  will  reveal. 

Table  1.  Partitions  of  Integers  1  through  5 


1 

2 

3 

4 

5 

1 

2 

3 

4 

5 

1+1 

2+1 

3+1 

4+1 

1+1+1 

2+2 

3+2 

2+1+1 

3+1+1 

l+l+l+l 

2+2+1 

l+l+l+l+l 
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5 


4+1 


3+2 


3+1+1  2+2+ 


l+l+l+l+l 


Figure  3.  Young  frames  for  partitions  of  5 


Figure  4.  Illegal  Young  frame 


Irreps  of  the  unitary  group 

The  group  U [n- 1)  is  contained  in  U  (n),  as  each  group  contains  the  generators 


of  the  smaller  groups.  As  an  example,  the  generators  of  U  (2) 


' 1  <P 

P)  P 

o 

P)  (P 

v0  0, 

v0  0/ 

V  oj 

5 

^0  P 

can  be  expressed  in  U (3)  as 


f  1 

0 

(P 

r° 

1 

(p 

ro 

0 

(P 

fO 

0 

(p 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

b 

10 

0 

b 

10 

0 

b 

10 

0 

b 

Thus  we  have  the  group  subduction  chain 


(121) 


(122) 
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C/(l)  cz  U  (2)  cz  C/  (3)  cz ...  cz  U  (n- 1)  cz  U  ( n ) 


(123) 


It  should  be  clear  that  the  unitary  group  U  (00)  contains  all  these,  (Jt/  («) . 

n 

Each  unitary  group  U (n)  is  a  continuous  group  and  will  have  uncountably 

infinite  elements.  Even  though  the  defining  representation  [34]  of  the  group  U  (n)  are 

the  n  x n  matrices,  it  can  be  represented  by  infinitely  large  matrices  with  an  infinite 
number  of  finite  irreps,  as  it  is  a  continuous  but  compact  group  [38],  The  irreps  of 

U(n)  can  be  labeled  with  Young  frames  of  at  most  n  rows  [38],  Figure  5  shows  the 

first  few  irreps  of  the  unitary  group  and  the  subgroups. 

The  GeTfand  Tsetlin  basis 
Weyl  tableaux 

In  the  UGA  we  will  use  Weyl  tableauxto  label  eigenfunctions  [38],  which  are 
Young  frames  whose  boxes  are  populated  with  tokens,  (usually  numbers,  but  it  must  be 
a  set  that  has  definite  ordering),  according  to  the  following  rules: 

1.  The  number  must  be  equal  to  or  greaterthan  the  number  to  the  left 

2.  The  number  must  be  strictly  greater  than  the  number  above 

In  the  UGA,  the  boxes  will  represent  electrons  and  the  numbers  will  represent  orbitals 
(or  spin  states),  thus  mimicking  the  indistinguishability  of  the  electrons.  Since  the 
numbers  represent  orbitals,  it  would  seem  this  system  allows  an  orbital  to  be  assigned 
to  more  than  two  electrons,  violating  the  exclusion  principle;  however,  it  will  be 
revealed  shortly  that  the  principle  of  antisymmetry  prohibits  such  Weyl  tableaux. 
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Figure  5.  Young  frame  representation  of  the  irreps  of  the  Unitary  Group 

Spin  States 

As  noted  previously,  the  spin  generators  are  those  of  U  (2) ;  from  Figure  5  it  is 

clear  that  the  spin  irreps  will  only  be  represented  by  one-  or  two-  row  frames.  Since 
there  are  only  two  spin  states,  spin-up  and  -down,  there  will  also  be  only  two  tokens,  1 
and  2,  or  a  and/?,  to  distribute  among  the  boxes.  Figure  6  shows  as  an  example  the 
possible  spin  states  for  a  four-  electron  system.  Note  that  each  box  in  the  top  row 
represents  adding  spin  1/2  to  the  total,  while  each  box  in  the  bottom  row  subtracts  spin 
1/2.  Thus  the  first  column  of  Figure  6,  the  totally  symmetric  irrep,  represents  the 
4x)<  =  2  spin  state  or  quintet;  the  second  column  represents  the  triplet,  and  the  third 
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5  =  2 


5=1 


5  =  0 


Figure  6.  Four-electron  spin  states 

represents  the  singlet.  Within  the  framework  of  the  total  spin,  the  system  can  have  a 
projection  onto  the  z-axis  from  +5  to  -5  in  integer  increments.  The  z  spin  component 
of  a  state  can  be  calculated  by  adding  1/2  for  each  1  and  subtracting  1/2  for  each  2,  and 
so  in  Figure  6  we  see  the  spin-structure  we  expect  out  of  a  four-electron  system.  Note 

that,  since  in  a  given  Weyl  tableau  for  U (2)  there  are  two  rows  and  two  numbers,  any 
column  with  two  rows  must  have  a  1  in  the  top  box  and  a  2  in  the  bottom.  Thus  the 
only  variability  in  the  ms  projection  is  within  the  single-box  columns.  For  this  reason, 
each  multiplet  does  not  have  a  unique  frame;  that  is,  any  number  of  two-boxed  columns 
could  be  added  to  the  Young  frame  and  the  same  multiplet  would  result.  See  Figure  7 
for  an  example  of  various  doublets.  If  we  were  only  interested  in  spin  states,  we  could 
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restrict  ourselves  to  single-row  irreps,  which  happen  also  to  be  the  irreps  of  SU  (2) 

[38];  nevertheless,  we  will  find  the  complete  U  (2)  irreps  invaluable  in  choosing  the 
correct  irreps  for  the  spatial  state,  as  we  will  see  shortly.  Furthermore,  we  can  use 
these  irreps  to  label  the  nodes  on  the  genealogical  construction  graph.  Figure  8  (see  also 
Appendix  I).  The  genealogical  construction  graph  shows  the  ways,  given  an  N  -  electron 
spin  function,  to  construct  N  +  1  -electron  spin  functions.  A  number  of  paths  (always 
stepping  to  the  right)  can  betaken  from  the  origin  to  each  node;  each  path  represents  a 
different  way  to  couple  angular  momenta  to  get  to  the  final  state.  While  the  Young 
frames  label  the  nodes,  they  give  no  indication  of  the  path  taken  to  that  node. 

Appendix  J  gives  the  relationship  between  paths  and  spin  functions. 


Figure  7.  Different  Young  frames  all  representing  doublets 


Spatial  States 

Recall  that  the  wavefunction  should  transform  as  a  member  of  an  irrep  of  a 
subgroup  of  U (2 n)  since  there  are  2 n  total  spin-orbitals.  Further,  if  we  restrict  the 

total  electronic  wavefunction  to  be  made  up  of  antisymmetric  tensor  products  of 
occupied  spin-orbitals,  it  must  reside  completely  in  the  totally  antisymmetric  irrep  of 

U (2 ft),  which  is  represented  by  a  single  column  of  N  boxes.  We  have,  however, 
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5/2 

2 

3/2 

spin 

1 

1/2 

0 

1  2  3  4  5 

electrons 

Figure  8.  Tableaux-labeled  genealogical  spin  construction  graph  [41] 

restricted  our  calculations  to  the  U ( n)®U (2)  subgroup  (by  restricting  spin  coupling), 
and  so  the  irreps  must  be  in  the  subduction  chain  from  the  antisymmetric  irrep  of 
U (2n)  to  that  group.  This  results  in  a  spin-dependent  Gel'fand-Tsetlin  basis.  This  will 
only  be  the  case  when  the  spatial  irrep  and  the  spin  irrep  are  mutually  conjugate  Young 
frames;  that  is,  the  frames  interchange  numbers  of  rows  and  columns  [38],  Figure  9 
shows  an  example  of  mutually  conjugate  Young  frames  that  would  result  in  a  legitimate 
antisymmetric  wavefunction.  Since  the  spin  irreps  can  have  no  more  than  two  rows, 
this  antisymmetry  condition  requires  that  the  spatial  irreps  can  have  no  more  than  two 
columns.  This  also  means  a  number  will  appear  no  more  than  twice  in  the  Weyl  tableau 
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spatial  irrep;  that  is,  an  orbital  will  be  assigned  to  no  more  than  two  electrons.  Since  the 


spatial  and  spin  Young  frames  must  be  mutually  conjugate,  the  spatial  Young  frame 
gives  the  spin  Young  frame,  and  so  it  will  be  unnecessary  to  display  both;  the  spatial 
frame  will  be  sufficient.  We  will,  however,  lose  the  z-projection  fo  the  spin,  which  will 
be  important  for  constructing  spin-orbit  wavefunctions  and  calculating  their  properties. 
Fortunately,  an  alternative  simple  accounting  system  for  ms  will  be  considered  later. 


Figure  9.  Mutually  conjugate  Young  frames 


An  example:  two  electrons  and  two  orbitals 

To  help  solidify  the  meaning  of  the  Weyl  tableaux,  a  simple  example  is  in  order. 
Consider  the  case  of  two  electrons  in  two  orbitals.  We  will  work  in  the  subgroup 
U (2)®U (2)  ciC/(4) .  The  possible  spin  states  are  a  singlet  and  a  triplet: 


or 

® 

where  we  have  retained  the  spin  frames  for  illustration.  For  the  singlet,  the  only  Weyl 
tableau  allowed  is 


a 

P 
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This  tableau  represents  a  totally  antisymmetric  spin  state  in  which  both  electrons  have 
been  given  different  spins.  This  is  the  familiar  function 

aP-Pa\  (124) 

where  position  indicates  electron  dependence  (i.e.,  the  first  function  in  a  product  is  of 
electron  1,  the  second  is  of  electron  2,  etc).  For  the  triplet,  there  are  three  Weyl 
tableaux: 

a  a 

These  tableaux  represent  totally  symmetric  spin  states.  The  first,  two  spin-ups;  the 
second,  a  mix;  and  the  third,  two  spin-downs.  These  correspond  to  the  spin  functions 

aa 

±[ap  +  pa\  (125) 

fifi 

Since  there  are  only  two  orbitals,  the  spatial  tableaux  are  the  same  as  the  spin  tableaux, 
and  represent  analogous  functions,  but  with  orbitals  rather  than  spin  states .  When 
paired,  they  result  in  the  following  states: 

Singlets  Triplet 

[ap-pa\±±  0  0  <^-^[<px<p2-(p2<px]aa 

0  0 

±[(pl(p2+(p2<P\][aP-pa]*]1\  2  |  0  0  |  a  |  P  |  <^\[(px(p2-(p2(p{}[ap  +  pa 

^<P2<Pi[aP-Pa\+*  0  0  ^^[cpxcp2-cp2cpx\pp 

0  0 
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These  are  precisely  the  antisymmetric  wavefunctions  we  expect  to  find  for  the  2- 
electron,  2-orbital  system;  thus  these  tableaux  label  a  spin-adapted  basis  for  the  system. 

In  this  example  the  antisymmetry  was  carried  completely  by  either  the  spatial  or 
the  spin  piece,  with  the  other  being  symmetric.  The  general  case  is  not  this  clean,  as 
pairs  of  tableaux  which  contain  neither  single  rows  nor  single  columns  are,  on  their  own, 
neither  symmetric  nor  antisymmetric,  yet  their  product  must  be  completely 
antisymmetric.  As  an  example,  a  3-electron,  2-orbital  system  has  two  doublets  (i.e.,  two 
paths  from  the  origin  in  Figure  8), 

f;j=[2  aa/3-(a/3a  +  /?««)] 


^[a/3a-  /3aa\ 

\$i[aPP-PaP] 


(126) 


corresponding  to  the  Weyl  tableaux 


a 

a 

and 

a 

(3 

P 

P 

none  of  which  have  definite  parity.  A  product  of  such  tableaux  is  not  as  simple  as  a 
product  of  a  spatial  and  a  spin  determinant,  but  rather  is  a  linear  combination  of 
determinants  which  is  antisymmetric.  Shavitt  has  developed  a  correspondence 
between  the  more  complicated  Weyl  tableaux  and  Slater  determinants;  [42]  may  be 
consulted  to  see  howto  express  CSFs  in  terms  of  Slater  determinants. 
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The  Subduction  Chain 


As  stated  previously  in  equation  (123),  each  unitary  group  contains  all  unitary 
subgroups  in  the  subgroup  chain.  Each  Weyl  tableau  represents  a  unique  subduction  of 

an  element  of  U  (n)  to  U  (l) .  The  subduction  of  a  tableau  in  U  (n)  to  U [n  -l) 
requires  removing  from  the  tableau  any  boxes  assigned  to  the  ntb  orbital.  Subsequent 
subduction  to  U (n  -2)would  then  remove  any  boxes  assigned  to  the  (n-\ )  "  orbital, 
and  so  forth  [38], 

Consider  the  following  Weyl  tableau  in  a  seven-electron,  seven-orbital  system. 


which  is  a  doublet.  The  tableau  resides  in  U  (7) ,  which  has  tableaux  of  no  more  than 
seven  rows.  To  subduce  to  U  (6),  we  must  remove  any  boxes  assigned  to  the  seventh 


orbital.  Since  there  are  none,  this  tableau  remains  unchanged  in  U[ 6).  Following  the 
prescription  above,  this  tableau  labels  the  subduction  chain  in  Figure  10. 

Paldus  Tableaux 

Each  frame  in  the  chain  can  be  described  by 

1.  the  number  of  two-element  rows  (  =  a  ) 

2.  the  number  of  one-element  rows  (  =  b  ) 

3.  the  number  of  zero-element  rows  (  =  c  ) 
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X 


1/(1) 


Figure  10.  Subduction  chain  for  an  example  doublet  Weyl  tableau  where  N=7 ,  n=7 


This  last  number  is  calculated  to  be  the  number  of  possible  allowed  rows  in  the 
subgroup  (i.e.,  n  rows  for  U ( n ) )  minus  the  number  or  rows  in  the  Young  frame.  Thus 

we  can  label  the  above  subduction  chain  with  the  numbers  a,  b,  and  c  in  order  of 
decreasing  orbital  number  (see  Table  2).  The  orbital  number  column  can  be  discarded 
since  it  is  the  sum  of  a,  b,  and  c,  and  we  create  the  more  simplified  table,  a  Paldus 
Tableau,  as  shown  in  Figure  11  [42],  The  Paldus  Tableau  is  equivalent  to  the  subduction 
chain  of  Figure  10,  and  is  hence  equivalent  to  the  original  Weyl  Tableau. 


Table  2.  Subduction  chain  from  Figure  10  in  table  form 


Orbital 

a 

b 

c 

7 

3 

1 

3 

6 

3 

1 

2 

5 

3 

0 

2 

4 

2 

0 

2 

3 

B 

B 

B 

2 

0 

1 

0 

0 

1 

0 

0 

0 

0 
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3  1  3 
3  1  2 
3  0  2 
2  0  2 
1  1  1 
0  1  1 
0  0  1 
0  0  0 

Figure  11.  Paldus  Tableau  corresponding  to  the  subduction  chain  in  Figure  10 

Gel'fand  Tableaux 

Finally,  we  can  alternatively  label  each  frame  in  the  subduction  chain  as  a  row  of 
numbers  whose  length  is  equal  ton  .  Starting  from  the  left,  we  insert  a  2  for  every  two- 
element  row,  followed  by  a  1  for  every  one-element  row,  ending  with  zeros.  This  is  very 
closely  related  to  the  Paldus  Tableau;  for  example,  the  top  row  of  the  Paldus  Tableau  in 
Figure  11,  indicating  three  2s,  one  1,  and  three  Os,  would  have  the  form  2  2  210  0  0. 

When  each  row  is  thus  labled,  the  result  is  a  Gel'fand  tableau.  The  Gel'fand 
tableau  in  Figure  12  is  equivalent  to  the  Paldus  tableau  in  Figure  11.  Thus  each  Weyl, 
Paldus,  or  Gel'fand  tableau  represents  a  specific  CSF.  A  Young  frame,  a  Paldus  tableau 
top  row,  or  a  Gel'fand  tableau  top  row  represents  a  family  (or  irrep)  of  CSFs  for  a  system 
of  n  orbitals,  N  electrons,  and  a  specific  multiplet. 

Ordering  tableaux 

All  of  thetableaux  in  this  irrep  can  be  ordered.  Each  tableau  set  has  its  specific 
ordering  rules,  but  those  for  the  Gel'fand  tableau  are  most  intuitive.  The  top  row  will  be 
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2  2  2  1  0  0  0 

2  2  2  1  0  0 
2  2  2  0  0 

2  2  0  0 
2  1  0 

1  0 
0 

Figure  12.  Gel'fand  tableau  corresponding  to  the  subduction  chain  in  Figure  10 

the  same  for  the  entire  irrep.  The  elements  of  each  subsequent  row  must  follow  the 
betweenness  condition;  that  is,  each  must  fall  inclusively  between  the  two  elements 
situated  in  the  row  directly  above  it.  The  Gel'fand  tableaux  can  be  ordered  according  to 
the  following  rule:  If  tableau  A  and  tableau  B  have  the  first  k  rows  in  coincidence,  then 

the  one  with  the  greater  number  of  2s  in  the  (£  +  l),/!  row  is  greater;  if  the  number  of  2s 
in  the  (k  +  lj1'  row  is  in  coincidence,  the  one  with  the  greater  number  of  Is  in  that  row 

is  the  greater  tableau.  Figure  13  shows  the  order  of  the  Gel'fand 
tableaux  for  the  system  of  three  orbitals  and  four  electrons  in  a  singlet,  followed  by  the 
Paldus  tableaux  and  Weyl  Tableaux,  including  the  subduction  chain.  They  are  labeled 
above  lexically.  The  basis,  thus  ordered  and  spin  adapted,  is  the  Gel'fand  Tsetlin  basis 
[38],  This  is  the  basis  COLUMBUS  uses  to  calculatethe  Cl  wavefunctions  and  the 
density  matrices  upon  which  analytic  gradients  and  DCTs  will  be  built,  and  shall  bethe 
basis  we  discuss  in  the  remainder  of  this  paper.  In  the  next  section,  we  will  discuss  how 
this  basis  is  used  in  the  graphical  approach  specifically  used  in  COLUMBUS. 
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Figure  13.  Gerfand  Tsetlin  Basis  for  W=4,  n=3,  singlet 


Distinct  Row  Tables 

Each  line  in  the  Paldus  tableau  format,  or  each  shape  in  the  Weyl  tableau  format, 
in  Figure  13  occurs  often  in  an  entire  l\l-electron  basis  where  n  »  N .  By  taking 
advantage  of  these  redundancies,  we  can  greatly  reduce  the  size  of  the  displayed  basis 
information  [43],  The  Distinct  Row  Table  (DRT)  contains  only  as  many  rows  as  there  are 
distinct  rows  in  the  set  of  Paldus  tableaux,  with  a  fixed  number  of  columns.  For  Figure 
13,  there  are  eight  distinct  rows;  while  there  is  not  much  savings  in  notation  (and  thus 
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disk  space)  at  this  level,  the  DRTs  of  larger  bases  of  realistic  problems  are  much  smaller 


than  their  Paldus  or  Weyl  tableaux  counterparts.  For  example,  a  30-orbital,  10-electron 
singlet  system  has  over  four  million  CSFs,  but  the  DRT  has  only  511  rows  [43],  Table  3 
shows  the  DRT  for  the  basis  in  Figure  13.  The  first  column  label,  j,  gives  the  index  of 
the  one-electron  basis  function,  which  is  also  the  unitary  subgroup  to  which  the  row 
belongs.  The  second  column  label,  k,  gives  a  unique  index  to  each  row.  The  next  three 

column  labels,  ak ,  bk ,  and  ck ,  label  the  row  itself  using  the  Paldus  step  numbers  (see 
Table  2). 

The  rest  of  the  DRT  essentially  shows  two  sets  of  data:  how  each  row  is 
connected  to  the  next,  and  how  each  row  is  lexically  ordered  with  respect  to  others.  To 
understand  the  first  point,  one  must  see  that,  in  addition  to  removing  an  orbital,  there 
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Table  3.  DRT  for  three-orbital,  four  electron  singlet  basis 
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are  four  possible  electron-related  actions,  indexed  by  the  step  number  d,  to  subduce 


one  row  to  another: 

1.  Remove  no  electrons,  spin  remains  the  same  (d- 0) 

2.  Remove  an  unpaired  electron,  spin  goes  down  (d- 1) 

3.  Remove  a  paired  electron,  spin  goes  up  (d-  2) 

4.  Remove  a  pair  of  electrons,  spin  remains  the  same  (d-  3) 

(Because  only  four  such  actions  are  available,  a  CSF  can,  in  addition  to  the  tableaux  in 
Figure  13,  be  identified  by  thetotal  spin  and  theorder  ofthese steps  taken  in  the 
subduction  chain.  This  will  be  useful  in  calculating  elements  of  spin-orbit  generators 
later  on.)  The  next  four  columns  of  the  DRT  show  to  which  row  the  given  row  subduces 
via  one  of  the  four  subduction  steps  listed  above.  These  indices  are  called  the 
downward-chaining  indices  [43],  For  example,  row  k  =  1  subduces  to  row  2  via  a  d  -  0 
step;  to  row  3  via  a  d  -  2  step;  and  to  row  4  via  a  d  =  3 step.  A  d  =  1  step  is  not  allowed 
from  row  1,  because  it  has  no  unpaired  electrons.  The  final  five  columns  are  labeled  by 
counting  indices.  The  final  column  of  the  DRT  gives  the  weight,  labeled  by  xk,  and 
counts  the  number  of  distinct  subtableauxthat  have  that  row  atthe  top.  The  y  indices, 
labeled  by  the  step  vector  d,  count  the  number  of  subductions  that  occur  lexically 
before  the  given  subduction.  There  is  no  column  for  y0k  since  the  d  -  0  step  always  has 
the  greatest  weight.  For  example,  row  1  subduces  to  row  2  via  a  d  -  0  step  in  one 
tableau.  Row  1  subduces  to  row  3  via  a  d  -  2step  in  two  distinct  Paldus  subtableaux; 
they  are  still  weighted  lower  than  only  one  tableau,  so  y2k  also  counts  1.  Finally  row  1 
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subduces  to  row4  via  a  d  -  3step  in  three  distinct  tableaux,  which  are  weighted  lower 


than  the  previous  three  tableaux,  and  so  receives  a  count  of  3.  This  can  be  seen  by 
comparing  the  DCT  in  Table  3  with  the  lexically  ordered  tableaux  in  Figure  13.  Note  that 
in  this  form  of  counting,  a  row  which  occurs  in  multiple  Paldus  tableaux  but  whose 
subduction  chain  is  identical  is  only  counted  once.  For  example,  the  subduction 
[l  0  0]  — >  [0  0  0]  occurs  in  three  different  Paldus  tableaux  in  Figure  13,  but  is  only 

counted  once  in  the  DRT. 

The  Shavitt  Graph 

While  the  highly-ordered  nature  of  the  DRT  is  useful  for  programming,  it  is 
difficult  to  comprehend.  The  Shavitt  graph  is  a  compact,  graphical  representation  of  the 
DRT  which  is  more  useful  for  human-readability  [43],  Figure  14  shows  the  Shavitt  graph 
corresponding  to  Table  3.  Each  row  is  represented  by  its  number  in  a  circle.  The  rows 


j  3 
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a=  2  11  0  0 


tail 


Figure  14.  Shavitt  graph  N= 4,  n= 3  singlet 
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k  are  distributed  in  a  grid  according  to  the  values  of  j,  a,  and  b.  To  the  left  of  each 
circle  is  the  x  counting  index.  Each  line  or  arc  represents  a  valid  subduction,  labeled  by 
its  arc  weight,  which  is  the  appropriate  y  index.  The  d  index  is  represented  by  the  slope 
of  the  arc.  Each  walk  from  tail  (bottom)  to  head  (top)  represents  one  of  the  CSFs.  Figure 
15  highlights  the  six  walks  that  represent  the  CSFs  for  the  example  basis.  Each  walk  is 
assigned  a  weight  by  the  sum  of  the  arc  weights  it  traverses;  thus  between  a  pair  of 
walks,  the  higher  weight  (lower  number)  is  assigned  to  the  walk  that  begins  with  the 
left-most  arc  [43], 

Representing  the  Non-Relativistic  Hamiltonian  Elements 

Elements  of  the  non-relativistic  second-quantized  Flamiltonian  are  calculated  as 

{<Pm  \  H  |  <Pn  )  =  Z  (Pm  I  K  I  <Pn  )  K  +  \  Z  I  \  <Pn  )  S rs,u  (127) 


Figure  15.  CSFs  of  N= 4,  n= 3  singlet  basis  represented  by  walks 
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where  the  bra  and  ket  are  CSFs  from  the  Gel'fand-Tsetlin  basis.  This  formulation  shows 


that  the  one-and  two-  electron  integrals  act  as  coefficients,  and  it  is  fundamentally  the 
representation  of  the  unitary  group  generators  in  the  Gel'fand-Tsetlin  basis  we  are 

calculating.  The  Etj  can  be  divided  into  three  types  of  operators: 

1.  Diagonal  or  weight  operators  where  /  =  j 

2.  Upper  triangular  or  raising  operators  where  i  >  j 

3.  Lower  triangular  or  lowering  operators  where  i  <  j 

The  CSFs  are  eigenvectors  of  the  weight  operators,  with  the  relationship 

<Pn)  =  #nmni  (128) 

where  nt  is  the  occupancy  of  the /"'orbital.  The  raising  and  lowering  operators  are 
adjoints  of  one  another,  so  we  will  only  consider  raising  operators  here.  There  are 
several  conditions  which  must  be  fulfilled  so  that  the  product  (cpm \Ejj  |y>n)  be  non-zero; 

the  equivalent  of  these  on  the  Shavitt  graph  is  that  the  two  walks  |  <pm)  and  |  q>\  must 
be  joined  from  tailtothe /-I  level  and  from  the  jth  level  to  the  head  [43],  Between 
these  two  levels,  the  walks  form  a  loop,  for  which  |^>H)must  always  be  to  the  right. 

Figure  16  shows  an  example  of  a  non-zero  loop  value  for  the  raising  operator  E2l  in  the 
example  basis.  The  walks)  2)  (solid)  and  1 3)  (hashed)  are  joined  at  levels  0  and  2  with 
1 2)  on  the  left  of  the  loop.  Their  common  walk  is  represented  by  the  dotted  line. 
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Figure  16.  Example  non-zero  unitary  generator  element 

To  calculate  the  loop  value,  the  loop  is  divided  into  segments.  Each  segment  is  a 
portion  of  the  loop  sitting  between  two  j  -indices,  e.g.,  [7-1,7],  [j,j  +  i\,...[i—l,i\.  A 
segment  is  given  an  index  k  which  is  equal  to  the  top  index  on  theShavitt  graph  or  the 
latter  in  the  pair  [  j,k].  In  the  example  in  Figure  16,  there  are  two  sections,  [0,l]  and 

[l,2],  shown  explicitly  in  Figure  17.  Figure  18  shows  the  seven  different  segment  types, 
Qk .  Just  like  the  generators,  the  segments  can  be  classified  as  weight  (W),  raising  (R),  or 

lowering  (L),  and  further  identified  as  being  unjoined  (no  overbar),  joined  at  the  top 
(overbar),  or  joined  at  the  bottom  (underbar).  Again  in  these  Figures,  the  bra  is  solid 
and  the  ket  is  hashed.  Additionally,  each  segment  can  be  identified  by  the  step  number 

of  the  bra  and  ket,  dk  and  dk,  respectively,  (where  the  primed  quantity  indicates  the 
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[U]  [0,1] 

k  =  2  k  =  1 


Figure  17.  Segments  of  the  loop  in  Figure  16 


Figure  18.  The  different  segment  types,  weight  (W),  raising  (R),  and  lowering  (L) 


bra),  and  the  values  of  b  at  the  top  of  each  segment,  bk  and  bk .  As  an  alternative  to 
this  last  requirement,  we  may  specify  the  difference  in  the  bk  values  in  addition  to  only 
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one  of  the  bk  values;  Shavitt  chooses  to  specify  bk  and  define  A bk  =bk  -  b\  [42], 
Segments  with  the  same  Qk,  d\,  dk,  and  A bk  have  congruent  shapes;  these 
parameters  are  collected  into  a  single  symbol 

Tt={Q,-,d'tdk-,Mt}  (129) 

called  the  segment  shape  symbol.  For  example,  the  shape  symbols  for  the  segments  in 
Figure  17  are 

Tx  = 

,  (130) 

T2  =  {/?;1,3;0} 

The  value  of  a  segment,  W (Tk,bk)  is  completely  determined  by  its  shape  symbol  and 

bk .  The  matrix  element  of  the  unitary  generator  between  the  two  CSFs  is  the  product 
of  these  segment  values: 

^\Et\m)  =  t[W(TlA)  (131) 

k=j 

The  values  of  W (Tk,bk)  are  given  in  Table  III  of  [42], 

Matrix  elements  of  generator  products  (eijkI )  are  calculated  in  a  similar  although 

more  complicated  manner  [42], 

The  spin-orbit  Hamiltonian  elements 

Adding  the  spin-orbit  operator  to  the  Flamiltonian  does  not  perfectly  fit  into  the 
above  calculations  for  two  reasons;  first,  the  generators  in  that  operator  are  not  spin- 
free;  and  second,  the  Shavitt  graphs  require  morethan  one  total  spin  irrep.  We  shall 
present  Yabushita's  solution  to  the  latter  problem  first  [11], 
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Multi-spin  Shavitt  graph 


When  introducing  spin-orbit  coupling,  our  example  system  is  incomplete; 
although  it  contains  singlet  CSFs,  these  could  now  couple  with  triplet  or  even  higher 
CSFs,  which  must  be  included  in  the  Shavitt  graph.  Figure  19  shows  this  graph,  with 
exclusively  triplet  arcs  and  nodes  marked  with  hashes  (see  Figure  14).  Unlike  the  single¬ 
spin  Shavitt  graph,  this  graph  now  has  two  heads,  one  for  each  possible  total  spin.  This 
means  that  FI amiltonian  elements  between  CSFs  of  different  spins  will  not  create  closed 
loops.  In  order  to  calculate  coupling  coefficients  between  CSFs  of  different  spins, 
Yabushita  devised  a  scheme  [11]  where  an  additional  index  is  added  to  the  top  of  the 
graph  and  the  closed-loop  formalism  is  used.  This  can  be  thought  of  as  adding  an 
additional  artificial  electron  and  orbital  to  the  Shavitt  Graph  (see  Figure  20),  creating  an 
irrep  of  U (n  +  l)that  will  subduce  to  both  the  singlet  and  thetriplet.  In  this  example, 
the  left  and  right  hashed  arcs  in  Figure  20  represent  the  subductions 


and 


respectively. 
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a=  2  21  11  000 

b=  1  02  10  210 


Figure  20.  Mixed-spin  Shavitt  graph  for  the  N= A,  n= 3  singlet  and  triplet  basis  with 

artificial  U(n+1)  head 
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Spin-orbit  generators 


Now  that  loops  between  CSFs  of  different  spins  are  built,  we  can  calculate  the 
matrix  elements  of  the  generators.  The  spin-orbit  generators  defined  in  (30)  are 
different  from  the  spin-free  generators  used  in  the  previous  subsection;  nevertheless, 
the  spin-orbit  generators  can  be  expressed  in  terms  of  the  spin-free  generators  [44],  To 
understand  this  process,  we  first  note  that  each  CSF  is  properly  identified  by  its 

subduction  steps  ( d  ),  total  spin  (s),  and  z-projection  of  spin  (ms).  As  mentioned 
previously,  the  subduction  chain  can  be  identified  by  the  order  of  steps,  d,  taken.  Thus 
we  can  identify  a  CSF  by  the  ket 


dSM 


(132) 


For  example,  the  CSF  |l)  in  Figure  13  may  be  represented  by 


|i}« 


\ 

d  = 

3 

13, 

S  =  0M  =  0  ) 

(133) 


(since  the  example  system  is  a  singlet,  the  z-projection  is  necessarily  zero.  For  non-zero 
spins,  the  physical,  two-column  Weyl  tableaux  do  not  specify  that  projection  without 
their  accompanying  two-row  spin  Weyl  tableaux).  Thus  the  matrix  element  of  the  spin- 
orbit  Flamiltonian  is  evaluated  as 


d'S'M' 


H 


d  S  m\  =  (d'S'M' 


r/usv 


dSM/K^y 


r,ju,s,v 


(see  equations  (19)  and  (31)).  Let  us  define 


(134) 
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(135) 


a  'ia 


(where  i  is  the  coordinate  of  the  i,h  electron),  so  that 


hso{i)  =  d(o)i)-q(i ) 


(136) 


(where  coj  is  the  spin  coordinate  of  the  i'h  electron),  even  though  in  the  actual 
implementation  the  matrix  elements  are  replaced  by  spin-orbit  potential  matrix 
elements  [45],  Since  we  know  that  q(i\  is  a  three-component  vector,  we  can  also  cast  it 


as  a  rank-one  spherical  tensor  by  defining  [11]  [46] 


4i(0  =  -^M0  +  M0) 
<7o  (0  =  clz  (0 

4i(0=\ffM0-M0) 


(137) 


(we  can  treat  cf.  equivalently),  whereby  the  spin-orbit  operator  can  be  recast  as 


A"(0  =  Z(-1)r«-r(0°'r(^) 

7 

Thus 

Klsv  ={r{i)/^{(oi)\hso{i)\s{i)v{(oi)) 

=  (r  (i )  //  (to. )  I X  (-1)7  q_r  (i )  crr  (». )  I  5  (/ )  v  ( ry,  )> 

=  <r 1  4-r  1  1  °V  1  V > 

r 

and  equation  (134)  becomes 


(138) 


(139) 
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d'S'M'  Hso  d  S  m)=  X  (d'S'M'  Er/uv  cl  S  M)£(-l)r  {r\cLy\ s)(ju\ cry  Ik) 


=  Z  (~lY  1  9-r  1  5>  (*' 'S  'M  '  \cjy\v)  cl  SM)  (140) 


{r\q_y\s)  (d'S'M  '  Zy{r,s)  dSM 


where  we  have  defined 


(141) 


which  is  itself  a  rank-one  spherical  tensor.  By  the  Wigner-Eckart  theorem  (see  Appendix 
I),  the  integration  of  that  tensor  can  be  cast  as  a  product  of  separate  geometric  and 
dynamic  pieces  [46]: 


r  s '  l  S  ^ 

d’S’M’\zr(r,s)\d  SM^  =  (-1)s"m'  ^  ^  (d'S'\\Z(i,j)\\ds)  (142) 


where 


S'  1  5 

-M'  y  M 


(143) 


is  the  Wigner  3-J  symbol  [46],  another  incarnation  of  the  Clebsch-Gordan  coefficients 


discussed  in  Appendix  I,  and  Id  '5 'll  Z(i,  j)  II  d  s\  is  known  as  the  reduced  matrix 


element  [46],  The  spin-orbit  Hamiltonian  can  be  expressed  as 


\S'-M’+r ,  ,  ..  ,  S'  1  S 


d'S'M'  HsodSM)='£(-lf  -  +r  {r  I  q_r  I  s)  (d'S'\\Z(i,j)\\d  S)  (144) 

r,s,r  \  JV1  y  JV1  J 


For  this  element  to  be  non-zero,  J  'must  be  obtained  from  d  by  substituting  orbital  r 
for  orbital  sin  J;  further,  y  =  M'-M  must  be  satisfied  [11];  thus  the  summation  is 


dropped: 
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d'S'M' 


H„ 


d  S M)  =  (-l)S  M  (H  qM_M.  I  s) 


S'  1  5 

y-M '  M'-M  Mj 


d'S'\\Z(i,j)\\ds)  (145) 


Drake  and  Schlesinger  [44]  show  that  the  reduced  matrix  element  (d'S'WZ  (/,  j)  II  dS 
can  be  expressed  as 

d'S'\\Z(i,j)\\ds)  = 

(-1) 


\$N+l+S  /  1 


[5’ 

5  1  I 

[  v, 

u 

1  s 

2  °N+l  J 

|  \d  N+l  ^N+ 1  MN+l 

^n+l,j^i,n+l  +  2  ^ ij 

d-N+ 1  ‘S'iV+l  -^iV+l 


(146) 


where  iV  +  1  is  the  index  of  the  extra  orbital,  n  +  l  is  the  index  of  the  extra  electron,  and 

fS'  S  1  ] 


1  ±  c 

2  2 


(147) 


is  the  Wigner  6-j  symbol  [46],  Thus  we  see  that  the  generators  are  now  spin-free  at  the 
expense  of  adding  an  orbital  and  an  electron  to  the  system. 

Now  we  proceed  to  simplify  the  matrix  elements  of  the  spin-orbit  Hamiltonian 
(equations  (144)  and  (146)).  Since  CSFs  can  couple  through  the  spin-orbit  operator  only 
if  their  total  spins  differ  by  0  or  ±1,  there  are  only  three  scenarios  for  SN+l : 

1.  it  is  a  half  integer  above  the  spins  of  equal-spin  CSFs, 

2.  it  is  a  half  integer  below  the  spins  of  equal-spin  CSFs 

3.  it  is  the  average  value  of  the  spins  of  CSFs  of  different  spins. 

For  example,  when  coupling  two  doublets,  SN+1  could  be  a  triplet  (scenario  1),  or  a 
singlet  (scenario  2);  when  coupling  a  singlet  to  a  triplet,  SN+l  must  be  a  doublet 
(scenario  3).  Given  these  three  scenarios,  the  coefficient  and  6-j  symbol  in  equation 
(146)  take  on  only  one  of  three  possible  values  [11]: 


70 


1. 


(S  +  1)(2S+1) 


2. 


S(2S  +  1) 
5  +  1 


(148) 


3. 


y/S+l 


(where  S  is  the  total  spin  of  the  ket)  respectively  pertaining  to  the  above  scenarios.  (In 
scenario  3  it  is  assumed  that  S '  =  S  + 1 ).  Yabushita  et  al.  [11]  define 


(f.) 

V  'J/N+ 1 


—  (d  N+l  SN+l  Mn+1 


(d  N+l  SN+1  Mn+ j 


En+ljEj  n+l  +  2  Ey 


SN+i  M  N+1 


p  +  J-  F 

c7,M+l,H+l,j  2  ij 


^iV+1  ^Ar+1  Mn+1 


(149) 


Just  as  in  the  previous  subsection,  these  matrix  elements  will  be  the  product  of  segment 
values.  The  quantity  can  then  be  expressed  as  the  product 


(150) 


The  segment  k  =  N  +  l  can  take  on  only  three  symbols,  pertaining  to  the  above  three 
scenarios: 


1. 

rM=!W;l,l;0} 

2. 

Vi  ={W/;2,2;0} 

(151) 

3. 

rM={RL;l,2;0} 

which  have  corresponding  segment  values  [11] 
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1. 


W(TN+l,bN+l)  = 


V  2(5  +  1) 


2. 


(152) 


3. 


W(TN+vbN+l)  =  1 


Thus  the  reduced  matrix  element  (d'S'\\Z{i,j)\\d  s\  can  be  simplified  in  each 


scenario  by  considering  the  product  of  ( F  )  with  the  appropriate  quantity  from 

'  v  y  Jn 

equations  (148)  and  from  equations  (152).  Considering  that  for  scenarios  1  and  2 
5  =  5 'and  for  scenario  3  5  +  1  =  5',  all  three  scenarios  result  in  the  same  product: 


d'S'\\Z(i,j)\\dS)  = 


5 '+5  +  1 


W, 


(153) 


In  turn,  the  elements  of  the  spin-orbit  Hamiltonian  are 


d'S'M' 


H 


dSM)  =  (-lf-M  (r  I  qM_M.  I  s) 


S'  1  5 

y-M '  M'-M  Mj 


5'+5  +  l 


(F«)„  <154> 


Note  that  all  references  to  the  additional  orbital  and  electron  have  been  removed  from 
this  form  of  the  spin-orbit  Hamiltonian  matrix.  This  method  of  calculating  spin-orbit 
generators  will  be  crucial  in  calculating  the  spin-orbit  DCTs. 

The  real  spin-orbit  Hamiltonian 

Molecular  systems  with  even  numbers  of  electrons  can  always  be  crafted  to 
result  in  a  real  Hamiltonian  [11].  Systems  of  odd  numbers  of  electrons  with  point 
symmetry  greater  than  D2  or  C2v  will  also  have  a  real  Hamiltonian  [47];  however,  odd- 


electron  systems  with  less  than  this  level  of  symmetry  will  have  a  necessarily  complex 
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Hamiltonian.  Yabushita  has  implemented  a  technique  which  adds  a  non-interacting 


electron  to  such  systems  so  the  Hamiltonian  will  be  real,  albeit  twice  the  dimension 
[11].  The  advantage  of  such  a  system  is  that  the  arithmetic  involved  in  diagonalizing  a 
real  Hamiltonian  is  less  computationally  costly  than  the  arithmetic  for  a  complex 
Hamiltonian,  despite  its  larger  size.  This  technique  also  fits  well  into  the  already- 
established  real  arithmetic  of  the  COLUMBUS  MRCI  programs  [7]  [8]  [9]  [10].  The 
disadvantage  of  the  extra-electron  technique  is  that  there  are  now  twice  as  many  states 
as  expected,  and  half  of  them  must  be  eliminated  as  they  do  not  represent  true  states. 
Kedziora  has  proposed  that  in  a  future  iteration  of  COLUMBUS,  the  extra-electron 
technique  be  abandoned  in  favor  of  a  smaller,  complex  Hamiltonian  with  the  correct 
number  of  eigenvalues  and  eigenvectors  [48], 
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III.  Formalism 


Normally,  energy  eigenvalues  and  wavefunctions  are  calculated  point-wise  in  the 
nuclear  manifold,  which  would  lend  calculation  of  gradients  and  DCTs  to  a  finite 
difference  method;  however,  finite  difference  methods  are  less  accurate  than  analytic 
methods,  especially  when  automated  [49],  Shepard  has  introduced  a  method  to 
evaluate  energy  gradients  of  Cl  wavefunctions  analytically  [3].  These  energy  gradients 
can  be  expressed  in  terms  of  atom-centered  basis  functions  and  geometry-dependent 
coefficients;  however,  due  to  Shepard's  method,  derivatives  need  only  be  calculated 
explicitly  for  the  original  atom-centered  basis  functions,  which  are  known  analytically. 
Lischka  et  al.  have  adapted  this  method  to  calculate  DCTs  analytically  [5],  which  is  the 
method  used  in  COLUMBUS.  We  shall  outlinethese  methods  here,  including  the 
modifications  necessary  to  accommodate  the  spin-orbit  Hamiltonian. 

Analytic  Spin-Orbit  Gradients 

Constructing  the  Solution  Wavefunctions 

Constructing  the  orthogonal,  geometry-dependent  basis 

Under  the  Born-Oppenheimer  approximation,  one  normally  chooses  a  point  R0 
in  the  nuclear  manifold  and  constructs  an  orthonormal  electronic  basis  suited  to  it.  At 
another  point  on  the  manifold,  R  =  R0+S ,  the  basis  at  R0  is  scrapped  and  a  new  basis 
is  created.  There  is  no  explicit  analytic  connection  between  these  two  bases.  Evaluation 
of  the  nuclear  gradient  on  electronic  functions  requires  knowledge  of  the  bases  at  R0 
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and  R  in  the  limit  that  8  — >  0 .  To  this  end,  it  behooves  us  to  construct  a  single 
geometry-dependent  basis,  rather  than  a  distinct  basis  for  every  point.  It  is  important 
to  keep  in  mind  that,  although  this  basis  will  be  constructed  to  apply  to  points  on  the 
manifold  within  a  S  neighborhood  of  R0,  because  of  the  limiting  process,  it  will  never 
actually  be  evaluated  at  any  point  except  R0 .  This  fact  will  aid  us  in  the  construction 
and  use  of  the  basis. 

In  the  following  derivation  by  Shepard  [3],  we  will  pass  through  many  bases.  We 
will  use  bracketed  superscripts  to  keep  track.  [  j]  will  indicate  the  atomic,  non- 

orthogonal  basis;  [C]  will  indicate  the  basis  that  is  geometry-dependent  and 
orthonormal  at  the  reference  geometry,  but  non-orthogonal  elsewhere;  and  [5]  will 

indicate  a  basis  which  is  geometry-dependent  and  everywhere  orthonormal.  Other 
bases  will  be  notated  as  they  are  introduced. 

For  a  given  geometry  R0 ,  we  have  a  set  of  normal  but  non-orthogonal  atom- 

centered  one-electron  functions  | Through  an  MCSCFstep  [35],  an 

orthonormal  set  of  molecular  orbitals  (RQ))  is  constructed,  with  the  relationship 

0(K)  =  C(X)Z(K)  (155) 

where 
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(156) 


(*.  Wl 

{^2(*o)| 

(*jW| 

Let  us  now  create  a  geometry-dependent  basis  in  a  similar  manner: 

^c|(K)  =  C(K,1)j(K)  (157) 

where  the  atomic  orbitals  are  allowed  to  vary  with  geometry,  but  the  coefficient  matrix 
is  not  reevaluated  at  every  point.  Such  a  basis  is,  however,  only  orthogonal  at  R0.  To 
remedy  this  shortcoming,  let  us  define  the  overlap  matrix  ,  such  that 

=  (158) 

In  the  [C]  basis, 

S[c]=C(R0)S[z]C\R0)  =  C(R0)z(R)®z\R)C\R0)  =  fic](R)<S>fic]t(R)  (159) 
Note  that  at  R0,  S^(/^,)  =  I.  At  every  other  geometry  R ,  there  exists  a 

transformation  ^(R)  such  that 

S[c]"^  ( R)S[C]  (i?)S[c]"^  (R)  =  I  (160) 

This  square  root  matrix  can  be  defined  by  the  expansion  [3] 

S[c]^(/?)  =  I-i(s[c](/?)-l)  +  |(s[c](/?)-l)2-...  (161) 

which  evaluates  to  I  at  R0 .  Also,  the  derivative  of  the  above  expression  evaluated  at 
the  reference  geometry  (which  will  soon  be  needed)  evaluates  to 
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dR  v  ’ 


*0 


_  At_i 

2 


Adcl 


S[CJ(tf) 


_d_. 

«o  ^  y 


+l(s[cl(«.)-i) 


-A, 

r„  dR  j 


(162) 


=  -\— s[c](^0) 

2  dR  V  07 


where  the  third  and  higher  terms  in  the  expansion  evaluate  to  zero  as  (Rq)  =  I 


With  this  transformation  matrix,  we  can  now  construct  the  basis 


^[s](7?)  =  SLcr/2(i?)^iCJ(i?) 


j[C]' 


(163) 


which  is  orthonormal  at  all  geometries  R  . 

Constructing  a  solution  away  from  the  reference  geometry 
We  have  created  an  orthonormal  one-electron  basis;  however,  the  actual  wave 
functions  will  be  represented  in  the  CSF  space  spanned  by  the  Gel'fand-Tsetlin  basis. 
Consider  the  following  highly  contrived  yet  didactic  example.  Suppose  we  are 
interested  in  a  system  of  six  electrons  with  four  orbitals  in  a  singlet  configuration.  Then 
there  are  10  CSFs  in  the  full-CI  space,  represented  by  the  Weyl  Tableaux  in  Figure  21. 
Recall  that  we  construct  the  Gel'fand-Tsetlin  basis  with  a  set  of  spatial  orbitals  from  a 


|1)  |2>  1 3)  1 4)  1 5)  1 6)  1 7)  1 8)  |9>  |10> 


Figure  21.  Gel'fand-Tsetlin  basis  for  N= 6,  n= 4  singlet 
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prior  MCSCF  step.  Suppose  in  the  MCSCF  step,  we  want  to  leave  the  first  and  second 


orbitals  doubly  occupied.  Then  we  are  restricted  to  the  MCCSF  space  of  |l(/?0)\(S> 
|2(/?0)^®  |5(i?0)^,  since  they  arethe  only  CSFs  in  which  this  occurs.  Assumethat  we 
already  have  the  MCSCF  solution  at  R0 .  Let  |mc(/?0)\  be  an  eigenfunction  of  H(i?0) . 
Figure  22  visualizes  this  eigenfunction  in  three-dimensional  CSF-space. 


|5(Ro)> 


Figure  22.  MCSCF  solution  vector  projection  in  3D  CSF-space 


Now  let  |^(f?))  be  an  eigenfunction  of  H (i?),  which  is  confined  to  the  space 

|l(i?)^(S)|2(i?)^(S)|5(/?)^ .  Since  both  of  these  functions  are  normal,  they  can  be 
transformed  one  into  the  other  by  a  unitary  transformation: 

|y/(tf))  =  exp(A(tf))|mc(/?0))  (164) 

where  A  is  the  antihermitian  product  of  parameters  and  unitary  generators 
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(165) 


MR)  =  Ta,AR)(E«-E,r) 

r<s 

These  two  wavefunctions  are  represented  as  coefficient  vectors  in  a  CSF  space.  They 
differ  in  two  important  ways:  first,  their  CSF  expansion  coefficients  are  different;  and 
second,  the  CSFs  themselves  differ  as  each  set  is  composed  of  MOs  that  have  resulted 
from  a  geometry-dependent  optimization  of  energy.  The  unitary  transformation 

exp(A(/?)jencompases  both  of  these  differences,  and  so  we  can  split  it  into  two 

separate  transformation  operators:  one  that  is  a  rotation  of  the  CSF  expansion 
coefficients,  and  one  that  is  a  rotation  (and  optimization)  of  the  CSFs  which  make  up  the 
basis, 

y/^  =  exp(^(/?)jexp(.P(/?)j|mc^  (^0))  (166) 

where  the  coefficient  rotation  exp acts  on  the  function  first,  followed  by  the  CSF 

rotation.  Flere  we  introduce  the  [/C]  basis,  which  is  geometry-dependent.  Figure  23 
shows  an  example  of  these  two  rotations. 

A  short  discussion  of  essential  and  redundant  orbital  optimization  is  now 
relevant.  The  MCSCF  step  rotates  the  orbitals,  which  consequently  rotate  the  CSFs.  This 
procedure  improves  the  solution  to  a  subsequent  (non-full)  Cl  calculation.  The  MCSCF 
rotations  can  be  partitioned  into  two  sets:  essential  and  redundant.  Essential  rotations 
are  those  that  result  in  a  lower  eigenvalue  in  the  MCSCF  step,  while  redundant  rotations 
have  no  effect  on  that  eigenvalue  after  a  reoptimization  of  the  MCSCF  CSF  coefficients. 
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|5(R„)> 


ep|mc(/?0)) 

|2(R„)> 


|i(Ro)> 


|5'(R)> 


|2'(R)> 


Figure  23.  (Left)  Rotation  of  a  vector  via  CSF  expansion  coefficients;  (Right)  Rotation  of 

a  vector  by  rotating  the  CSFs 


Redundant  rotations  take  place  in  invariant  orbital  subspaces.  The  subsequent  Cl  step 
will  also  rotate  CSFs,  but  only  via  coefficients— not  the  orbitals.  It  too  will  have  a 
partition  of  essential  and  redundant  rotations.  If  the  MCSCF  invariant  space  is 
contained  in  the  Cl  invariant  space,  any  set  of  MCSCF  invariant  rotations  will  have  no 
effect  upon  the  Cl  eigenvalue;  however,  if  some  component  of  the  MCSCF  invariant 
space  lies  outside  the  Cl  invariant  space,  seemingly  arbitrary  rotations  in  the  MCSCF  step 
will  produce  different  energies  in  the  Cl  step.  For  this  reason,  an  additional 

transformation,  exp(z(/?)j ,  known  as  orbital  resolution  is  required  to  ensure  the  Cl 

energies  are  well-defined,  and  thus  there  is  a  well-defined  limiting  process  to  obtain  the 
derivatives.  Let  us  continue  the  example  to  see  how  this  transformation  affects  the 
solution  wavefunction. 

Since  we  have  left  orbitals  1  and  2  doubly  occupied  in  this  example,  the  orbital 
subspace  spanned  by  those  orbitals  is  invariant.  (For  further  discussion  of  invariant 
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subspaces,  see  [3].)  We  can,  therefore,  arbitrarily  rotate  orbitals  1  and  2  by  the  unitary 


operator  exp(z(/?))  in  the  invariant  subspace  and  the  CSFs  will  still  span  the  same  CSF 

space  (see  Figure  24).  That  is,  a  solution  vector  in  the  basis  on  the  left  side  of  Figure  24 
can  be  expressed  in  the  basis  on  the  right  by  a  rotation  of  the  CSF  expansion 

coefficients,  exp  ^-Z(i?)j,  that  exactly  offsets  the  rotation  of  the  CSFs  in  the  invariant 

space  (see  Figure  25).  Thus  we  have 

y/[z]  (R'fj  =  exp (Z (tf))exp(-Z (R)J |  y/[K]  {Fifj 

=  exp^(7?)jexp^Z(f?)j^exp^-Z(f?)jexp^(f?)jJ|/?7c^^  (^0))  (1^7) 

=  exp  (^))exp(z(f?))exp(p(f?))|7«c[s](7?0)^ 

where  exp(-Z(i?))  as  been  rolled  into  exp^P(i?)j,  as  they  are  both  rotations  of 
coefficients  only.  The  two  vectors,  while  physically  the  same,  have  distinct 


|5'00> 


Figure  24.  (Left)  the  MCSCF  optimal  CSFs;  (Right)  An  arbitrary  rotation  in  the  invariant 
space  leads  to  an  equivalent  optimal  set  of  CSFs,  where  the  energy  is  the  same  as  the 

previously  optimized  energy 
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|5'(K)> 


Figure  25.  Different  representations  of  the  same  vector  due  to  a  redundant  rotation 


representations.  For  this  reason,  we  have  introduced  the  [z],  or  resolved  basis,  as  our 
final  one. 

Now  suppose  that,  for  our  Cl  step,  we  relax  the  double  occupancy  restriction  of 
orbital  2,  but  we  keep  orbital  1  doubly  occupied  and  require  at  least  one  electron  to  be 
in 

orbital  3  at  all  times.  The  Cl  space  is  |r(/?))®|2'(/?)^®|3'(/?)^®|6'(/?))®|8'(/?)^  . 

The  Cl  solution  will  not  contain  the  MCSCF  solution  at  R  because  |5'(/?)^  is  absent,  but 


it  will  contain  the  projection  onto  the  2 '(R))  space,  which  we  explicitly 


highlight: 


|C/(«))  =  [exp(z(«))|r(«))<l'(«)|exp(-Z(«))]|C/(«)) 
‘exp(z(R))|2'(R)>{2'(R)|exp(-Z(R))]|C;(R)) 

£  M*))W*)|c/(*)) 


(168) 


mgMCSCF 
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Unless  ^5'(i?)|exp(z(i?))|C/ (i?)^  =  0,  the  eigenvalue  will  not  be  minimal,  and  will 

depend  on  the  operator  K  chosen.  Thus,  while  Z  is  redundant  in  the  MCSCF  step,  it 
must  be  chosen  such  that  ^5'(/?)|exp(z(/?))|C7 (R)j  =  0  in  orderto  minimize  the  Cl 
energy. 

Given  the  correct  orbital  resolution,  the  Cl  solution  vector  at  R  can  be  expressed 
not  only  as  a  rotation  of  the  MCSCF  solution  at  R0 ,  but  as  a  rotation  of  the  Cl  solution 
vector  at  R0  as  well.  Thus  the  Cl  solution  at  R  can  be  expressed  as 

^c/](i?))  =  exp(^(i?))exp(^(i?))exp(JP(i?))|^]  (R0))  (169) 

All  of  the  nuclear  geometry  dependence  has  now  been  removed  from  the  wave  function 
and  resides  in  the  optimization  operators. 

The  Analytic  Gradient 

We  can  now  use  the  Cl  solution  wavefunctions  at  an  arbitrary  nuclear  geometry 
to  construct  the  analytic  gradient.  The  Cl  solution  vectors,  which  are  vectors  of  CSF 
expansion  coefficients,  satisfy  the  (electronic)  matrix  Schrodinger  equation 

H(R)ci(R)  =  s1(R)ci(R)  (170) 

where  s1  is  the  eigenvalue  of  the  Cl  wavefunction  y/, ,  represented  here  by  a  vector  of 
coefficients  in  the  CSF  space,  c1 .  In  this  case,  H  is  the  electronic  Flamiltonian  matrix 
which  includes  both  non-relativistic  and  fine-structure  contibutions,  including  the  spin- 
orbit  Flamiltonian  (see  equation  (12)): 

H(R)  =  H„(R)  +  H,„(R)  (171) 
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We  can  take  the  derivative  of  the  Schrodinger  equation  and  left -multiply  it  by  cJ : 


i'  (*)£(«(*)?'  («)) = c1  ( R )1'  ( R ))  • 


(172) 


Applyingthe  chain  ruleto  both  sides. 


as, 


(£')c'(*)c'(*)  +  £'(*)cy(*)^(e'(*)) 


(173) 


The  first  term  on  the  right-hand  side  is  the  delta  function  c 7  ( R )c7  ( R )  =  SM ,  yielding 


_d_ 

dR, 


(174) 


For  the  gradient,  we  will  assume  that  /  =  J  ;  later  we  will  makethe  alternate 
assumption  for  the  DCTs.  We  can  now  isolate  the  energy  gradient. 


_d_ 

SR, 


{s’  (*))  =  ?'  (*)i(H(/i))c'  (R)+s’c’  (R)£(c’  (R)j-s’c’  (R)^(c’  (R)) 
=  S’(R)i(H(R))cJ(R) 


(175) 


Let  us  use  the  following  shorthand  notation  for  equation  (175): 


£*(k)  =  W*)IhWM*)) 


(176) 


where  the  superscript  a:  denotes  differentiation  with  respect  to  a  nuclear  coordinate. 

Although  this  equation  appears  to  be  the  Helmann-Feynman  Theorem,  it  is  not. 
The  Helmann-Feynman  Theorem  applies  to  exact  wavefunctions,  whereas  this  formula 
applies  to  Cl  wavefunctions,  which  are  only  exact  within  the  space  spanned  by  the  basis. 
Following  equation  (171),  let  us  split  equation  (176)  into  two  pieces: 


£'(R)  =  (^(R)|H„(fl)>,(R))  +  <^(7;)|H„(fi)>,(R)) 


(177) 
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Shepard's  formalism  already  addresses  the  first  term  on  the  right-hand  side  of  equation 
(177)  [3]  (also  see  Appendix  G  for  more  detail),  and  so  we  will  focus  on  the  last  term, 
following  in  a  parallel  manner. 

Per  Shepard's  formalism,  the  spin-orbit  gradient  in  the  resolved  [z]  basis  can  be 
expressed  in  the  orthogonal  [5]  basis  as 

(yf]  WK'W  K1  w) = (iT1  to)K'  W  |rJs)(«o)) 

+  (rfl(^)|[H^l(S„),Z~(S„)-]|^1(j;a))  (178) 

+  (^!Sl  («„ )  |  [H™1  («o )  •  K'“"  («„ )'  ]  |  («„ )) 

where  higher  orders  of  commutators  are  considered  negligible,  and  the  bra-ket  notation 
indicates  vector-matrix  multiplication  (cf.  equation  (269)  in  [3]).  This  result  stems  from 
the  form  of  the  solution  wavefunction  at  an  arbitrary  geometry  as  expressed  in  equation 
(169).  Let  us  first  address  the  first  term  on  the  right-hand  side  of  equation  (178),  then 
subsequently  we  shall  address  the  final  two  terms,  which  are  treated  similarly  to  one 
another. 

The  First  Term:  The  Orthogonal  Basis  Derivative  Term 

Using  the  methods  of  second  quantization,  the  first  term  on  the  right-hand  side 
of  equation  (178)  can  be  reduced  to 

(rf  («»)|h“(7!0)>1s|(R.))  =  Ec!?’K1(^)|  e;.~  I V?  (K )) 

r/iSV  (179) 

=  7>  (h*|!|“  («<j!|  (fi„ )  |  %rm  |  (T1  («,,))) 
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where  E  is  a  non-spin-averaged  unitary  group  generator,  and  the  matrix  hso 
represents  the  operator  hso  defined  to  be  proportional  to  the  spin-orbit  operator: 


l,m 


im\ 


(180) 


(cf.  equation  (11);  see  [11]).  The  derivative  of  the  spin-orbit  integral  matrix  is  not 
calculable  by  the  present  form  of  any  integral  program  in  the  COLUMBUS  suite; 
however,  these  terms  are  available  in  NWCHEM  in  C,  symmetry  [50],  When  this 
operator  is  integrated  over  electronic  coordinates  between  spin-orbitals,  we  see  that 


h[s]so  =  /r[s] 


rjusv 


5[s]v\  =  (r^ 


s  v 


E  {r[s]  |X^/  (rJ  (4 )  I lm)  (lm  4[s] )  {m  I 

X  l,rn 


sxlV 


(181) 


Let  us  define  a  new  non-spin-averaged  density  matrix : 


ErMSV  ^JK))  =  EcK  W5J(r:/?o)  Erflsv  {r :  RQ))  (182) 


V  /  i 


(cf.  definition  (34))  where  are  the  CSFs  and  c\  are  the  corresponding  wavefunction 
coefficients.  Combining  equations  (182),  (181),  and  (179),  we  can  rewrite  the 
orthogonal  basis  derivative  term: 


(4%)  H \so\R0)x  ysf]{R0))  = 

££('%  :^)|£[6  W(4)l  H(lm\s[s](re  :Ro))'Z(M\Sz\v)'Lcti  {^]{re:Ro)\ErMsy\<plS]{re:R0) 


(183) 


which  we  will  separate  into  a  spin-independent  integral  piece: 
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and  a  spin-dependent  density  piece: 


11V  ij 


which  we  will  call  the  new  spin-orbit  density  matrix.  Thus  we  see  that  the  first  term  on 
the  right-hand  side  of  equation  (178)  reduces  to  the  trace  of  the  matrices 


(T1  K)|H“  («„)’  |tT'  (*„))  =  7>(,M  (*,,)'  ■  ZM  (*„)) 


(186) 


As  a  result  of  this  dissertation,  the  spin-orbit  density  matrices  are  now  available  in  the 
COLUMBUS  program  CIDEN,  and  are  traced  with  the  derivatives  of  the  one-electron 
integral  matrices,  qz,  in  our  modified  version  of  NWCHEM.  (In  NWCHEM,the  integral 
matrices  and  their  accompanying  analytic  gradients  are  multiplied  by  /  to  force  them  to 
be  real.) 

The  Second  and  Third  Terms 

The  nuclear  dependence  will  now  be  dropped  for  brevity.  In  the  second  and 
third  terms  on  the  right-hand  side  of  equation  (178) ,  all  three  matrices  have  a  second- 
quantized  form: 

u[s]  =  y  hso[s]E 

so  r/isv  r/usv 

rjusv 

Z“=2>:(4-£„)  (lg7) 

rs 

k . =2>r(4-£„) 


where 


<J?=ZZWr>|s>> 


(188) 


leading  to  the  forms 


rP|[HS,.Zf"*]H’l)‘SC*(»f,|  E  ''?£l[^>.v(E„-E„)]|riS|)  (189) 


rlil|[Hl;l,r“]|rfi)=XC'(rl'1|  E  *isi-[^>.v'(^.-4)]|^sl)  (wo) 


Wemust  now  determinethe  meaning  of  the  commutator  Er,fts,v,^Ers  -  Esr\  ,  which 


appears  in  both  terms.  Using  the  commutation  rule  for  unitary  group  generators, 

\E  E  1  =  E  8vp  -  E  8pa 

rpsv  ’  tpua  rfj.ua  st  svtp  ru 


(191) 


we  conclude  that 


/«v  /  /v  ^  \  _  n  n 

E.  .  ,\E  -E  )  =E.  Sva-E.  8pa  +  E ,  RSvp -E ,  „SMp 

r  fis  v  ’  \  rs  sr  I  r  fisa  s  r  s  vra  r  s  r  fisp  s  r  s  vrp  r  s 

—  E  Sva  +  E  8pa  -  E  8vp  +  E  8mP 

^ r'  prau s'  s  ^  s'vsau r'r  ^  r' fir  pu  s' s  .9  'vsfiur'r 


(192) 


When  the  requisite  sum  over  spin-orbitals  is  applied  and  the  spin-orbit  integrals 
included,  we  find 

V  hso  ,  \e  .  ,  ,(e  -E  )]  =  y/iso  E.  -Yhso,E, 

r  fis  v  r  fis  v  ’  \  rs  sr  J  r  fira  r  fisa  sas  v  s  vra 

r'jus'v  r'  fi  s'v 

+  YhS0  E  - y hso  E  -Yhs°  E  (193) 

^  Z-i  r’firft- ^r'LisP  sfls'v ^s'vr/3  /_,nr' psa^r' pm  I 

r'  fj  s'v'  r' n 

+  ^jKas'v^s'vsa  ~^jK'fis/3^r'nr/3 

s'v  r'  fi  sV 

Since  h50is  symmetric,  we  take  advantage  of  the  fact  that 


1  SO  1  so 

Ksv  =  KrM 


(194) 
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as  well  as  the  freedom  to  choose  independent  summation  indices  to  rewrite  equation 


(193) as 


Xhso  Te  (e~E )]  =  Y(  hso  E  +hso  E  )  -  V  ( hso  E  +  hso  E  ) 

lr' jus'v  ^r'fis'v’y^rs  «' /  J  £j\  ’’rat/i  tftsa  ^  "rflt^tfis/3  J  "salfl ^tpra  T  ls/3tfl ^tfir/3  ) 

r'jus'v  t/i  tji 

Hhso  E  +hso  E  )  +  'Y(hso  E  +hso  E  ) 

ri,sj3tju1^tjurj3  ’^sat/u^t/ira  J  'at ft1"' t jusa  r  Lr  (ft  ju1^  t  jus  (!  J 


(195) 


which  further  reduces  to 


V  hso  ,  \e,  .  ,(e  -E  ]]  =  2Y  hm,  E,  -  IT  hm,  E.  (196) 

r  jus  v  r  /us  v>\  rs  sr  J  rvtju  tjusv  svt/u  tjurv  v  • 


Each  of  theterms  on  the  right-hand  side  of  this  equation  is  an  element  of  a  spin- 
contracted  Fock  matrix  for  the  spin-orbit  Hamiltonian.  (See  Appendix  G  for  definitions 
of  non-relativistic  Fock  matrices).  Let  us  define 


(^slL=Z'-:.»K1|4.>isl 


=ZZ(4sl),,(4sl)„ 


(197) 


(/-). yK1) -(4?), 


(198) 


so  that  equations  (189)  and  (190)  yield 


'!s|  [H'i'.z-"] 


Zmcx  /  r  so  \  _  -+mcx  ^rsc 

^ rs  orb  J  ^  J orl 


(199) 


rf 1  |[h1'1.  k""]|  r!11} = XC  (fZ ),,  =  *“  •  fZ 


(200) 
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respectively,  yielding 


(T1  (K )  I H™  (R„  Y  |  vf]  K  )>  =  Tr  (qW  (R, )  ■  zf  (R„ ))  +  z""  •  /"  +  k' ”  •  /”  (201) 


Orbital  Resolution  Vector 

It  has  been  shown  by  Shepard  that  the  form  of  the  orbital  resolution  vector,  z  , 
depends  on  the  method  of  resolution  [3];  if  several  methods  are  used,  the  contribution 
from  each  method  must  be  summed.  The  three  methods  available  in  the  MCSCF 
program  are  natural  orbital  (NO)  resolution,  F-Fock  matrix  resolution,  and  Q-Fock  matrix 
resolution;  however,  the  F-Fock  matrix  resolution  is  not  available  for  use  in  gradients 
and  DCTs  at  this  time.  The  other  two  methods  result  in  the  following  forms  of  the 
orbital  resolution  vector: 

1.  Natural  Orbital 


~NOx 
** rs 


Z) 


[K]x 


d[k]-d[k] 


(202) 


2.  Q-Fock  Matrix 


;Qx 


(203) 


where  the  natural  orbital  resolution  uses  the  one-electron  non-relativistic  transition 
density  matrix  defined  in  equation  (35),  and  the  Q-Fock  resolution  uses  the  Fock  matrix 

elf1 -2^42(2*2 -£1)DM  (204) 
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The  Product  zx  ■  f 

While  the  product  of  the  orbital  resolution  vector  gradient  with  the  orbital 
gradient  vector  will  remain  the  same  for  the  non-relativistic  term,  the  spin-orbit  term 
will  require  a  few  modifications  based  on  the  results  presented  above.  That  product  can 
be  expressed  as: 


■/i=2Ed1 


d[k]  -  d[k] 


=2X]'K)„ 


(205) 


7„*  =  2£et 


e'f'-e'f1 


=Ze!f'“(4«)r 


(206) 


(The  undefined  diagonal  elements  are  resolved  by  defining  A  to  be  zero  on  the 
diagonal.)  Note  that  these  new  A*  matrices  differ  from  Shepard's  formalism  only  in 
that  they  have  the  spin-orbit  Fock  matrices  in  the  numerator,  while  all  other  matrices 
remain  the  same.  The  result  is  that  these  terms  are  evaluated  parallel  to  Shepard's 
formalism,  with  the  exception  that  the  A*  matrices  are  substituted  for  the  non- 
relativistic  variety. 

The  derivative  of  the  Q-Fock  matrix  (definition  (204))  is 

ejf1' = zMf1' + 1(2*2*  -  *2*K' + Z(**2  -  *2 )  w) 

tu  tu 

When  partially  transformed  into  the  [S]  basis,  this  derivative  is 


2/,»  +  2{/1W,r}ri  +  E(2£'-£l>W 

tu 

+E(2{*|SU1 -{*[JU*} . )d«+E(2  «2-*2)o2* 

tu  rS  U  r  SU  tu 


(208) 
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where  {A,B  }  =  AB  +  B  A .  Three  types  of  terms  exist  in  this  gradient: 

1.  Terms  involving  derivatives  of  integral  matrices 

2.  Terms  involving  derivatives  of  K 

3.  Terms  involving  derivatives  of  transition  density  matrices. 


The  natural  orbital  resolution  will  only  involve  the  third  type  of  term  because  it  has  no 
Fock  matrices. 


Q-Matrix  Terms,  Type  1 

These  terms  have  the  form 


(2°-  ■  /„*)' =2>Lsl‘K  L + Z(2«2.'  -  l 

rs  rstu 

=  7>  (h[SI*D£  ) + }  7>  (g[sl'd° ) 


(209) 


where  we  have  defined 


-i(AS)„oW 


Q-Matrix  Terms,  Type  2 

These  terms  have  the  form 


(210) 


(2-'7i)"=E2{A|sU-l 


=  rr({/.M,K'}Di)+irr({«|!1,K'}d«) 

=  rr(K-’  (h^Dj  +  Djh[sl  +  ig|s|d»  +  ig|s|d»  ))  (211) 

-  -27V(K'Ff ) 


=  k 


mcx 


-pQso 
J  orb 
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where  we  have  analogously  defined 


F,;-hWDj+ig|s|dJ  (212) 

Q-Matrix  Terms,  Type  3 

These  terms  have  the  form 

( ?e'v  •  fork  )  =  Z ( 2S[mu  -  sill )  D^]x  ( A°  )ra 

rstu 

tu 

=IErf(»P1(aOJ*.+*«)hw)  (2i3) 

n  tu 

=EK(»ri  |2Qihw> 

n 

=px-?S° 

where  (n^  is  the  complement  space  of  the  MCSCF  solution  vector  rnc^j ,  and  we 
have  defined 

(oi)„  sZ(2*-  -*2,)K)„  (2i4) 

tu 

Natural  Orbital  Terms,  Type  3 

These  terms  have  the  form 

(z”‘  •/,:)"'= ZfW-  (a;  i 

rs 

n  rs 

=  TjPn{n±]\2A°o\mc[K]) 

n 

=  y .  fDso 

—  r  n  J  csf 

Combining  equations  (215),  (213),  (211),  and  (209)  into  equation  (201),  we  have  the 
more  refined  form  of  the  spin-orbit  gradient: 


(215) 
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wf] I | vT1}  =  7>(hw'D»  )  +  \Tr (gM'< )  +  Tr (qf  •  Z« ) 


(^Xt  pVt) 


fSO  _|_  ?fiw 

J  orb  J  orb 

~rQso  ~rDso 
Jcsf  J  csf 


(216) 


Matrix  Turnover 


We  can  simplify  the  final  term  in  equation  (216)  by  defining 


Xmcx=(kx  t  px'\ 


/'total  _ 

so 


^cso  ,  ~rQso 
J  orb  J  orb 

/-*'Qso  ,  rDso 
csf  Jcsf 


(217) 


Xmcx  is  known  as  the  first-order  response  of  the  MCSCF  wave  function,  and  is  shown  by 
Shepard  [3]  to  be  equal  to 


j^mcx  Q  J  orb  /  \  ^ 

dRx  7*  me  V  / 
<J  csf  J 


(218) 


at  the  reference  geometry,  where  Gmc  is  the  Hessian  matrix,  fff  is  the  orbital  gradient 
vector  and  fff  is  the  CSF  gradient  vector.  Thus  we  have  that  the  last  term  in  equation 


(216)  is  equal  to 


fSO  .  fQso 
J  orb  J  orb 


Jorb  J  orb  =  frcxQmc-l  J,o,al 
V  *  /  -fQso  |  rDso  J  sc 

>  J  csf  J  csf  J 


*  csf  J  csf 


(219) 


To  the  right-most  matrix-vector  product  we  assign  the  symbol 


y^so  _  ^^mc—\j?U 


(220) 


which  we  separate  into  orbital  and  CSF  pieces: 


n  mexfy  me— 1  rtotal  rmex  j  so  .  rmex  n  so 

^  ^  J  so  ~  Jorb  *  \rb  +  Jcsf  *  \sf 


(221) 
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Analogously  to  Shepard's  method,  we  assert  that  the  orbital  piece  has  the  form 


J?  ■ =' Tr  (h|!)’D«  )  +  h  Tr  (gw’d» ) 

where 


(222) 


(223) 


and  A“feis  the  matrix  form  of  the  vector  As0°b.  Similarly,  we  have 

fT -K;  =rr(h^)  +  irr(g^dt)  (224) 

where 


"  (225) 

(C  )rsm  =  )n  1  in±  I  Ksru  +  Csnu  +  Ksut  +  e.rut  \  mC) 

n 

Substituting  equations  (224)  and  (222)  into  (216),  we  find  that  the  spin-orbit  gradient  is 

<=rr(h[sl'(Dj+D»+Di))  +  irr(gW'(d“+d*+di))  +  rr(qW’-zW)  (226) 
which  we  simplify  as 

<  =  rr(h[s1'  (d:  ))  +  irr(gp1'  (d“ ))  +  rr(q|sl*  -Z[s| )  (227) 

This  form  of  the  spin-orbit  energy  gradient  must  be  transformed  back  to  the 
atomic  orbital  (AO)  basis  so  that  it  can  be  traced  with  the  analytic  AO  basis  gradient 
integrals.  Note  that  the  corrections  due  to  orbital  resolution  result  in  traces  with  the 
non-relativistic  integral  matrices  rather  than  the  spin-orbit  integral  matrices.  This  is  due 
to  the  fact  that  the  spin-orbit  integral  matrices  were  not  involved  in  the  MCSCF  step. 
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For  this  reason,  D“f  and  d'“  can  be  summed  up  with  the  non-relativistic  effective 
density  matrices  and  traced  accordingly.  This  fact  makes  coding  the  spin -orbit  effective 
density  matrices  in  the  program  CIGRD  as  simple  as  adding  the  spin -orbit  Fock  matrix  to 
the  non-relativistic  Fock  matrix  before  any  further  calculations  are  performed. 

Transformation  to  the  Atomic  Basis 

The  three  traces  in  equation  (227)  are  in  the  [5]  basis,  and  must  be  transformed 

to  the  atomic  basis.  The  density  matrices  do  not  change;  however,  the  integral  matrices 
transform  as: 


h[s]x 


(228) 


(where  the  braces  indicate  the  anticommutator)  leading  to  the  transformation  of 
equation  (227): 


Shepard  has  also  shown  that  the  fourth  and  fifth  trace  operations  in  equation  (229)  are 
equal  to 


^7r({h[c],S[c]'} 

jyo?  [c]  \ 

SO  J 

\  =  Tr\ 

^[c]-*fLsc[c]  j 

i7>({g[c],S[c]*' 

K[c|) 

\  =  Tr\ 

Js[C].vF2.vC[C]j 

(230) 


where  and 


jji2sc[C] 


are  effective  Fock  matrices  defined  as 
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(f»'|c|LsZ^K|c|)b 

t 

(F,r|c|)  -z««(<cw)  U3D 

tuv 

(f,:|c|)  -(Cc|)  h(f,tic]) 

\  /  rs  \  '  rs  \  /  rs 

In  a  similarfashion,  the  last  trace  operation  in  equation  (229)  becomes 

\  Tr  (|q[c],  SA[C] }  •  Z[c]  j  =  Tr  (sf[c]Fs[oc] )  (232) 

This  Fock  matrix  was  defined  in  equation  (197).  Substituting  equations  (232)  and  (230) 
into  (229),  the  spin-orbit  gradient  is  defined  completely  in  terms  of  traces  of  products  of 
integral  and  density  matrices  in  the  [C]  basis,  whose  form  is  equivalent  in  the  atomic 
basis.  After  the  transformation  of  the  matrices  from  the  MCSCF  molecular  orbital  basis 
to  the  AO  basis,  the  spin-orbit  energy  gradient  becomes 

e  xso=  Tr(hx[z]  (vZ[x]))  +  \Tr(g[x]  (<C  W))  +  7r(q*W  -Z1 W)-7>  (s*W  (Cw  +F^))  (233) 
These  terms  are  added  to  the  non-spin-orbit  form  of  the  gradient, 

<  =  Tr  (hx[z]DM'W  )  +  i  7r(gxWd'"'W  )  -  Tr  (s*WF  ^  )  (234) 

(see  AppendixG)  to  yield 


In  the  actual  implementation,  the  standard  Fock  matrices  F  and  Fvo  are  added  together 
early  in  the  program  CIGRD.  The  effect  of  this  combination  is  to  eliminate  the 
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distinction  between  the  various  spin-orbit  Fock  and  density  matrices  defined  in  this 


section  and  their  spin-contracted  counterparts;  that  is, 

j ~ytOt  ^  j -ytOt  _ ^  j -ytOt 

dtof+d'°;  ~^dlor 

17  .  ~w?tot  .  -|7 

F  tot  '  r  so  '  rtot 

thus  simplifying  the  definition  of  the  energy  gradient: 


(236) 


(237) 


Analytic  Derivative  Coupling  Terms 

Elements  of  the  derivative  coupling  term  matrix  P  have  the  form 

/*(*)' =(M*)|4tM*)>  (238) 

(where  we  are  using  Lischka's  notation  for  the  DCT  [5]).  The  wavefunctions  i//,  are 
expanded  in  the  CSF  space  as 

Mr:-R))=Zci'Wk(r:-R))  (239) 

i 

Using  this  form,  the  elements  of  the  DCT  matrix  are 

/,.(*)’  =  Z4  (*)(«>,  ('"’)!«('- : «)))  <240) 

ij 

Applying  the  chain  rule, 

fn  W  =  Ik  (R)(<P,  k  :  (*))|  W  (r  ■ «)) 

(241) 
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In  the  first  term,  the  derivative  operation  on  the  Cl  coefficients  does  not  participate  in 


the  integration,  leaving  only  the  delta  function  (cp.  (r :  R)\(p,  (r :  R'j)  =  cT  .  In  the  second 

term,  the  Cl  coefficients  again  exit  the  integration,  but  the  gradient  operation  remains. 
This  leaves 

/»W’=Sc/W4r(c-'(/i))+Sc/'(*WW(pJ('-:J')l4rk.('-:*)>  <242> 

i  ij 

The  first  term  is  a  vector  dot-product,  and  emphasizes  the  change  in  the  Cl  coefficients; 
thus  it  is  called  by  Lischka  et  al.  the  Cl  DCT 

/“  (RY  =  (R)  =  cJ  (R) £c‘  (R)  (243) 

i 

The  second  term  emphasizes  the  change  in  the  CSFs,  and  is  called  the  CSF  DCT: 

f,7r(RY  (*W  R)\ikv  ( R ))  (244> 

ij 

The  addition  of  the  spin-orbit  operator  to  the  Flamiltonian  will  significantly  alter  the  Cl 
DCT,  which  will  be  discussed  below.  In  COLUMBUS,  the  CSFs  are  formed  from  the 
orbitals  defined  in  the  MCSCF  step  which  does  not  account  for  the  spin-orbit  operator. 
For  this  reason,  the  CSF  DCT  is  unchanged  by  the  addition  of  the  spin-orbit  operator  at 
the  Cl  step.  Lischka's  formalism  for  the  CSF  DCT  is  reviewed  in  Appendix  H . 

The  Cl  DCT 

Recall  that  we  found  the  energy  gradient  by  assuming  I  -  J  in  equation  (174). 

To  solve  for  the  Cl  DCTs,  we  make  the  alternate  assumption,  i.e,  I  ^  J  .  Applying  this 
assumption  to  equation  (174)  yields 
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(*)&(<?'  (*))  =  (*'  (*)"«'  (*))"' Z1  Wi5(H(S))f'  (R)  (245) 

As  we  did  with  the  energy  gradient,  we  can  split  the  Hamiltonian  into  two  pieces. 


fS  («Y  =(*-(«)-*'  («))"  (*)  £(H.  (R))c‘  (. r ) 

+  («'  («)-4J  (*))■'  (*)4r(H„  (R))c'  (R) 


(246) 


and  define 


/“”(«)'  -  c1  (*)4-(H0  (*))?'  (*) 


(247) 


Thus  the  Cl  DCT  is  very  similar  to  the  Cl  energy  gradient,  with  the  notable  exception  that 
it  has  an  energy  difference  in  the  denominator,  and  the  coefficient  vectors  are  from 
different  (rather  than  the  same)  wavefunctions.  For  this  reason,  the  formalism  for  the 
analytic  spin-orbit  Cl  DCT  exactly  mirrors  that  of  the  analytic  spin-orbit  energy  gradient, 
with  the  addition  of  the  energy  difference  in  the  denominator  and  the  substitution  of 
transition  density  matrices  in  place  of  the  standard  matrices: 


Ae  fr  =  Tr  (hx[z]  (d“/,[z]  ))  +  \  Tr  [gx[x]  (d ™JI[z] )  j 

+  Tr  (qxW  •  ZJI[z] )  -  Tr  (F*0cy/[zI  +  F'/W 


(248) 


This  quantity  is  added  to  the  standard  Cl  DCT  (see  equation  (407)  in  Appendix  G)  and  the 
standard  CSF  DCT  (see  equation  (432)  in  Appendix  H)  to  produce  the  full  form  of  the 
DCT: 
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Ae  /“*CSF*”-'  =  Tr  (h’VI  j  D"w  +  D';;"'1'1  +  D”"'W ))  +  7r(q*w  ■  Z"w) 

+  ±7>(g'W  (d,;W  +dFf  ■ "M  +d:J'M))  (249) 

-rr(s*1(F*1+F"'W  +  F,-l'w))+A4rr(D"W‘(;i0)/°FW(;i0)') 

In  the  actual  implementation,  just  as  with  their  standard  counterparts  in  the 
formulation  of  the  energy  gradient,  the  spin-orbit  density  and  Fock  matrices  are 
combined  with  their  spin-contracted  counterparts, 


ny/M  +  T\totJIU] 

tot  so 


^DJ,[z] 

tot 


Fj/M  +  vscJI[x]  v  YJI^ 

A  tot  '  A  so  T  A  tot 


(250) 


to  yield  a  cleaner  formulation  of  the  DCT: 

As  f«+CSF+s°*  =  Tr (hA[z]  +  D^F//[z] ) j  +  Tr  (q't[z]  •  Z//[z] ) 

+{rr(g'W(d^4d™w))-r,'(s-w(F';b4F:w))  (251) 

+  ASrr(D"W"(R0)/°fU(4)'j 
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IV.  Results  and  Discussion 


Test  Case  and  Implementation 

In  order  to  validate  the  formalism  presented  in  Chapter  III,  we  have 
implemented  it  in  a  hybrid  of  NWCHEM  and  COLUMBUS  codes,  and  have  chosen  K  He 
with  a  Stuttgart  basis  [51]  [52]as  a  test  case.  K  He  was  chosen  was  because  of  its 
relative  simplicity  for  speedy  calculations,  the  availability  of  a  spin-orbit  potential, 
evidence  of  an  avoided  crossing  (where  DCTs  should  be  significant),  and  its  possible 
application  to  a  Diode-Pumped  Alkali  Laser  (DPAL).  Energy  surfaces  for  the  excited 
states  of  K  He  are  shown  below  in  Figure  26. 
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Figure  26.  K  He  MRCI  energy  surfaces 
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Calculation  of  spin-orbit  gradients  and  DCTs  requires  derivatives  of  spin-orbit 


integrals  and  potentials.  COLUMBUS  typically  uses  one  of  two  integral  programs: 
ARGOS,  which  is  capable  of  handling  spin-orbit  potentials  and  integrals;  or  DALTON, 
which  can  produce  gradients  of  atomic  integrals.  Since  neither  program  currenty  does 
both,  NWCHEM  has  been  leveraged  to  produce  the  spin-orbit  integrals  to  be  fed  into 
COLUMBUS,  as  well  to  calculate  thespin-orbit  integral  gradients  used  in  thetraces  in 
equations  (235)  and  (249).  Kedziora  has  modified  NWCHEM  to  writethe  integrals  into 
the  Standard  Integral  File  System  (SIFS)  format  used  by  all  COLUMBUS  programs.  The 
disadvantage  of  using  NWCHEM  for  integrals  is  that,  unlike  DALTON  or  ARGOS, 
symmetry-adapted  integrals  are  not  available.  This  shortcoming  requires  that  all 
calculations  must  be  done  in  the  Ci  symmetry  group,  thus  increasing  calculation  time, 
memory,  and  space  required. 

The  current  COLUMBUS  code  that  produces  the  MRCI  wavefunctions,  including 
MCDRT,  MCSCF,  MCUFT,  MOFMT,  CIDRT,  TRAN,  CISRT,  and  CIUDG,  has  not  been 
modified.  The  density  program,  CIDEN,  has  been  modified  to  correctly  interpret  the 
multi-headed  spin-orbit  DRT  or  Shavitt  graph  (see  Figure  20)  and  to  produce  the  anti¬ 
symmetric  spin-orbit  density  matrices  Z  in  equation  (185).  CIGRD,  which  produces  the 
effective  density  matrices  and  the  Fock  matrix,  has  been  modified  to  create  the  spin- 
orbit  Fock  matrix  Fm  defined  in  equation  (197).  This  matrix  is  then  added  to  the 
standard  Fock  matrix  defined  in  equation  (365)  of  Appendix  G.  This  sum  effectively 
combines  the  matrices  defined  in  equations  (205)  and  (206)  into  those  defined  in 


103 


equation  (380),  as  well  as  combining  F“  with  Frof  in  equation  (235).  A  separate 
symmetric  version  of  Fm  is  also  supplied,  as  required  by  equation  (235).The  spin-orbit 
density  matrices,  which  do  not  require  an  effective  counterpart,  are  passed  through  to 
the  CIGRD  output  files.  After  transformation  to  the  atomic  basis,  the  symmetric  density 
and  Fock  matrices  are  passed  back  into  NWCHEM,  ratherthan  DALTON,  to  evaluate  the 
atomic  integral  gradients  and  perform  the  final  traces.  For  DCTs,  the  antisymmetric 
(non-spin-orbit)  density  matrix  is  still  passed  to  DALTON  to  evaluate  the  final  trace 
operation  in  equation  (249). 

The  Density  Matrices 

Modification  of  the  density  matrices  for  use  with  spin-orbit  DRTs,  as  well  as  the 
creation  of  the  spin-orbit  density  matrices,  is  the  cornerstone  of  the  implementation. 

We  can  validate  these  matrices  by  using  them  to  calculate  the  energy  expectation  value 
of  the  wavefunctions: 

(y/i\H\y/I)=Tr(hB)  +  ^Tr(gd)  +  YJTr(qzZz)  (252) 

X 

(see  equations  (32),  (184),  and  (185)).  Notethat  this  is  also  the  trace  of  the  total  Fock 
matrix: 
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r''(F)=I  Frr 

p 

= xf  S- KD; ^EEK),  (z, ), 

p  V  '  tuv  x  1 

=  Z  V  +  2  Z  ^  Amv  +  Z  Z  K  )„  ( : ZZ  )„ 

pt  ptuv  X  Pl 

=  rr(hD)  +  irr(gd)  +  X?'--(q/Zi) 


(253) 


This  value  should  be  nearly  equal  to  the  energy  eigenvalues  calculated  by  CIUDG,  which 
approximately  solve  the  Schrodinger  equation: 


(254) 


where  is  the  approximate  eigenvector.  The  difference  between  these  values  is 

dependent  upon  the  convergence  tolerance  of  the  MRCI  eigenvalues. 

Figure  27  compares  the  MRCI  eigenvalue  surfaces  with  the  Fock  energy  surfaces 
for  the  test  case  for  the  p-manifold  excited  states.  Point-wise  differences  are  on  the 
order  of  1(T10  hartree  or  smaller.  Each  of  these  states  is  normally  doubly  degenerate, 
but  due  to  the  addition  of  the  ghost  electron  [11],  each  state  is  quadruply  degenerate. 

It  is  sufficient  to  specify  the  energy  or  gradient  of  only  the  lowest  of  each  quadruplet . 

Since  the  effective  density  matrices  differ  from  the  standard  density  matrices  in 
that  they  account  for  rotations  in  the  MCSCF  invariant  space,  the  energy  expectation 
value  they  construct  should  also  be  approximately  equal  to  the  MRCI  eigenvalues. 

Figure  28  compares  these  values  with  the  MRCI  eigenvalues.  Their  difference  should 
now  depend  not  only  on  the  convergence  of  the  MRCI  wavefunctions,  but  also  upon  the 
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Figure  28.  K  He  Effective  Fock  -  MRCI  energy  difference 
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convergence  tolerance  of  the  MCSCF  Hessian  inverse.  Point-wise  differences  are  in 
general  on  the  order  of  10~10,  but  grow  an  order  of  magnitude  when  approaching  the 
repulsion  wall.  Since  energies  themselves  are  only  converged  to  10  6  hartree  in  this 
calculation,  these  small  energy  differences  are  most  likely  numerical  artifacts  due  to  the 
different  methods  used  to  calculate  either  energy. 

Analytic  Gradients 

Analytic  gradients  have  merit  on  their  own,  in  that  they  can  be  used  for 
geometry  optimization  algorithms  [53],  Adding  this  open-shell  spin-orbit  capability  to 
analytic  gradients  has  the  potential  to  improve  such  optimizations.  While  this  is  a 
contribution  in  its  own  right,  our  interest  in  them  is  that  the  code  required  to  produce 
the  gradients  is  a  very  large  part  of  the  code  that  produces  the  DCTs.  The  successful 
calculation  of  these  gradients  partially  validates  the  calculation  of  the  DCTs. 

The  linear  geometry  of  K  He  allows  for  easy  calculation  of  a  finite-difference 
gradient  to  which  the  analytic  gradient  can  be  compared.  Figure  29  through  Figure  31 
compare  finite-central-difference  gradients  using  a  0.2  A  split  to  analytic  gradients. 
Differences,  for  the  most  part,  are  less  than  1%. 

Derivative  Coupling  Terms  and  the  Adiabatic  Mixing  Angle 

Effects  of  the  Arbitrary  Rotation  of  Wavefunctions 

DCTs  suffer  an  additional  complexity  that  the  gradients  do  not.  Each  MRCI 
wavefunction  may  have  an  arbitrary  phase  factor,  which  does  not  affect  the  energy 
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Figure  29.  Finite  vs.  analytic  gradient  of  K  He  II 17  state 
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Figure  30.  Finite  vs.  analytic  gradient  of  K  He  Tlv  state 
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Figure  31.  Finite  vs.  analytic  gradient  of  K  He  state 
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expectation  value  or  its  gradient.  Let 

y,)  =  e-‘°\i7,)  (255) 

Then  for  an  operator  A 

{y'i \a\y,)  =  eW(v, \Ae~w\¥l) 

=  ewe-wtyI\Ae~u>\y/I)  (256) 

=  (y, \A\Vt) 

and  so  the  phase  factor  is  superfluous.  If  we  instead  consider  a  transition  property 
between  different  wavefunctions,  we  find  that  the  phase  does  indeed  matter: 
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(v,\a\ Vj)  =  eie'  (Vi  \Ae~Wj  \y/ ,) 

=  ew,e~i9j  (i//,  \Ae~w\^j)  (257) 

=  e,{6,~0j)(i//I  \A\y/j) 

Since  COLUMBUS  only  uses  real  wavefunctions,  the  arbitrary  phase  is  limited  to  ±1,  but 
it  will  still  require  careful  deliberation,  especially  as  the  DCT  crosses  the  axis. 

In  addition  to  the  arbitrary  phase  angle,  there  will  also  be  an  arbitrary  mixing 
angle  between  degenerate  wavefunctions;  the  doubling  of  degeneracy  due  to  the  ghost 
electron  only  compounds  this  issue,  creating  a  quadruple-degeneracy.  To  understand 
how  this  mixing  angle  affects  the  DCTs,  consider  an  arbitrary  vector  in  the  space 

spanned  by  the  E17  states: 

/2 

V)  (258) 


We  can  project  this  vector  into  the  space  of  gradients  of  the  Ylv  states: 

/2 


I 

j= i 


nj/Vnf 


_8_ 

dR 


^>=Ziftni)z(n; 


7=1 


i= 1 


K.)(A 


v 


From  this  equation  we  conclude  that 


i= 1 


(259) 


(260) 


that  is,  the  DCTs  are  simply  the  projections  of  the  gradients  of  the  Tly  states  onto  the 

E17  space.  Let  us  apply  an  arbitrary  unitary  transformation  to  the  basis  of  the  space 

/2  /2 

(for  all  practical  purposes,  we  treat  the  space  and  its  dual  equivalently): 
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_e_ 

dR 


/=! 


(261) 


While  the  gradient  wavefunction  itself  has  not  changed,  the  DCTs  have  been  changed  by 
the  transformation.  This  means  that,  unfortunately,  a  DCT  such  as  /nT 


Zy),  being  a 


projection  onto  a  single  arbitrary  basis  function,  is  not  well-defined.  The  Euclidean 
norm  of  the  gradient  wavefunction,  however,  is  well  defined.  That  is, 


Thus  we  see  that  the  physically  meaningful  quantity  is  the  norm  of  the  DCTs.  This  norm, 
as  we  shall  presently  show,  is  the  true  DCT  for  a  rotated  basis. 


From  the  four 


n,,  wavefunctions  to  the  four 


wavefunctions  there  are  a 


total  of  16  different  DCTs,  represented  within  the  matrix 
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(263) 


2/y 


In  reality,  pis  an  off-diagonal  submatrix  of  the  entire  DCT  matrix  for  all  wavefunctions; 


however,  if  we  are  only  interested  in  the  coupling  between  the  jnT  [  space  and  the 
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E^j  space,  we  need  only  consider  the  matrix  formed  between  these  eight 

wavefunctions  and  their  duals,  for  which  the  diagonal  blocks  are  zero: 

(  0 

P  =  t 

l-P  °J 

If  we  assume  that  p  is  not  singular,  then  there  exists  a  similarity  transformation  by 
which  it  may  be  diagonalized: 


prf  =  v  pv 


(264) 


(265) 


f  0  prf] 

V-1 

r 

r  o 

V-Prft  0, 

V  0 

V-Pf  °y 

v  0 

It  follows  that  the  off-diagonal  blocks  of  the  matrix  P  can  be  diagonalized  by  the 
transformation 

P'  = 

We  now  find  P'  in  a  new  basis  in  which  the  original,  arbitrarily  rotated  and  phased 
wavefunctions  have  been  re-rotated  to  yield  a  canonical  set  of  wavefunctions  {\u'v_ 

and 


(266) 


E^M.  This  transformation  has  reduced  the  number  of  DCTs  from  sixteen  to  four 
(the  four  eigenvalues  of  the  off-diagonal  block),  each  of  which  couples  only  one  Yl'y 

wavefunction  to  only  one  E^  wavefunction.  Considering  that  now  each  row  and 
column  of  P'  contains  only  one  non-zero  element,  equation  (260)  tells  us  that  each 

E',/\,  and  that  the  length  of  the  derivative 


m  )  is  collinear  with  a  distinct 


wavefunction  is  the  DCT.  Since  the  wavefunctions  in  each  respective  space  are,  for  our 
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purposes,  indistinguishable,  there  should  be  no  distinction  of  the  coupling  between  any 
two  pairs,  and  we  should  expect  all  four  new  DCTs  to  have  the  same  magnitude. 

As  an  example  of  the  preceeding  discussion,  consider  the  DCTs  of  the  K  He 
system  at  5.0  A  (9.4  bohr).  The  calculated  DCTs  are  found  in  the  matrix 


6.68  x  ICT6 

-2.42  x  10~6 

-1.72  xlO-4 

5.11x10 

-5.10  x  1CT2 

1.63  xlO"5 

-6.51  xlO-6 

5.33x10 

-5.54  x  10~6 

1.04  x  10  3 

5.11  xlO~2 

1.75x10 

-1.84  xlO"5 

—5.11 x 10  2 

1.04  x  10~3 

1.56x10 

which  is  diagonalized  by  the  matrix 


^3.47x10  3  -4.09x10  \  4.80x10^+4.08x10  \ 

3.47  x  r3  +  4.09  x  10~'  i  4.80  x  10_1  -  4.08  x  10_1  i 


V  = 


7.06  x  10 


7.06  x  10 


2T 

6 

-4 

6 

J 


-3  -1  '\ 

9.88x10  '  -4.09x10  i 

9.88  x  10"3  —  4.09  x  10  1  /' 


(267) 
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2.89x10  1  +5.00x10  'i  1.02x10  3  -7.34x10  \  -2.89x10  '  +5.00x10  */ 
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-1 
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-1 
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(268) 


When  the  matrix  P  built  from  the  matrix  in  equiation  (267)  is  subjected  to  the 
transformation  in  equation  (266),  the  resultant  DCT  matrix  is 


0  5.11  x  10  -  5.05  x  10  /  0  0 

0  0  5.11  x  10  '  +5.05  x  lo"4  i  0 

0  0 


0  -2.55  x  10  -  4.43  x  10 


0 


0 


0 


-5.11x10  -5.05x10 


0  -5.11x10  +5.05x10  i  0 


0  2.55  x  10  -  4.43  x  10 


0  2.55  x  10  +  4.43  x  10 


(269) 


While  each  of  the  DCTs  in  P'  is  complex,  we  have  the  option  of  effecting  an  arbitrary 
phase  rotation  between  each  pair  of  wavefunctions,  allowing  us  to  remove  the 
imaginary  portion  of  each  DCT: 
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0  0  0  0  -0.0511  0  0  0 

0  0  0  0  0  -0.0511  0  0 

0  0  0  0  0  0  -0.0511  0 

ooooo  o  o  -0.05H  (270) 

0.0511  0  0  0  0  0  0  0 

0  0.0511  0  0  0  0  0  0 

0  0  0.0511  0  0  0  0  0 

0  0  0  0.0511  0  0  0  0 

As  hypothesized,  we  see  that,  in  the  canonical  basis,  the  DCTs  between  each  distinct 
pair  of  wavefunctions  are  indeed  of  the  same  magnitude.  From  this  exercise  we  also 
conclude  that  we  need  only  calculate  four  DCTs,  not  sixteen.  From  the  four  distinct 
elements  of  any  given  row  or  column  of  the  DCT  matrix,  we  can  calculate  a  Euclidean 
norm,  which  is  the  DCT  in  the  canonical  basis. The  radial  DCT  in  the  canonical  basis, 
calculated  as  the  norm  of  the  DCTs  in  the  first  row  of  P,  is  shown  in  Figure  32. 

Adiabatic  Mixing  Angle 

This  exercise  has  shown  that,  despite  the  eight  wavefunctions  involved,  we  are 
effectively  seeking  the  adiabatic  mixing  angle  for  only  two  wavefunctions  and  so  we  can 
calculate  the  adiabatic  mixing  angle  according  to  equation  (103).  Figure  33  shows  the 
result  of  that  integration  (as  integrated  from  R  =  oo).  Figure  34  shows  the  non-adiabatic 
surfaces  and  the  resultant  mixed  diabatic  surfaces.  Figure  35  shows  the  off-diagonal 
diabatic  coupling  surface.  The  strength  of  the  coupling  surface  depends  upon  both  the 
size  of  the  coupling  angle  and  the  energy  difference  between  the  two  adiabatic  surfaces; 
thus  the  coupling  surface  continues  to  grow  well  past  the  area  where  DCTs  are  strong, 
into  a  region  which,  for  most  interactions,  would  be  energetically  disallowed.  These  are 
the  surfaces  that  may  be  used  in  nuclear  dynamics  calculations. 
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Figure  34.  Diabatic  vs  adiabatic  K  He  surfaces 


Figure  35.  Diabatic  K  He  coupling  surface 
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Validation 


We  have  two  sources  of  validation  of  our  DCT  calculations.  First,  since  the 
analytic  gradient  algorithm  is  such  a  large  part  of  the  DCT  algorithm,  the  correct 
gradients  in  the  previous  section  speak  to  the  fidelity  of  the  DCTs  calculated  with  the 
same  code.  Second,  Werner  established  a  method  to  estimate  the  coupling  angle  using 
the  MRCI  wavefunction  coefficients  [12]: 


^  Z 

where  the  Euclidean  norm  is  taken  of  the  coefficients  of  the  eigenvector,  ci  72 , 

which  coincide  with  CSFsthat  have  II 17  symmetry  (for  simplicity,  only  the  reference 

/2 

CSFs  are  considered  in  an  MRCI  calculation).  We  can  generalizethis  method  to  spin- 
orbit  wavefunctions  by  using  the  Euclidean  norm  of  all  coefficients  pertaining  to  CSFs 
whose  physical  occupation  is  appropriate.  Figure  36  shows  our  calculated  mixing  angle 
compared  to  the  angle  estimated  via  the  Werner  method.  Similarly,  Figure  37  compares 
our  DCT  calculation  with  the  derivative  of  Werner's  angle.  The  Werner  angle  matches 
the  angle  calculated  with  the  DCTs  by  less  than  0.02  radians,  an  acceptable  error  given 
the  rudimentary  nature  of  the  approximation. 
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Figure  36.  Coupling  angle  as  calculated  via  DCTs  vs.  calculated  via  Cl  coefficients 
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Figure  37.  DCT  as  calculated  vs.  derivative  of  Werner's  angle 
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V.  Conclusion 


A  tool  has  been  developed  to  calculate  energy  gradients  and  DCTs  of  spin -orbit 
wavefunctions  at  the  MRCI  level.  While  several  approximation  methods  have  been  in 
place,  no  method  has  been  available  that  actually  calculates  these  quantities  for  spin- 
orbit  Cl  wavefunctions  until  now.  This  work  makes  such  a  method  available,  having 
been  aggregated  to  the  MRCI,  DCT,  and  spin-orbit  methods  already  in  place  in 
COLUMBUS. 

With  this  formalism  in  place,  spin-orbit  energy  gradients  and  DCTs  can  now  be 
analytically  calculated  for  large-atom  systems  where  spin-orbit  contributions  to  the 
Hamiltonian  are  non-negligible,  such  as  open-shell  systems  like  the  Alkali-noble-gas 
mixtures  that  may  be  used  in  a  DPAL.  While  the  energy  gradients  will  be  useful  in 
geometry  optimization  problems,  the  DCTs  will  provide  diabatic  potential  energy 
surfaces  for  nuclear  dynamics  calculations. 

KHe 

We  have  calculated  such  surfaces  specifically  for  the  K  He  system.  The  maximum 
derivative  coupling  for  K  He  appears  to  occur  well  outside  the  energy  well  of  the  Yly 

manifold  (13.2  bohr  separation),  indicating  that,  despite  its  small  magnitude  (0.117 
radian/bohr),  the  derivative  coupling  may  play  a  non-negligible  role  in  the  dynamics, 
and  should  be  taken  into  account  for  a  relevant  DPAL-related  analysis.  Additionally,  the 
coupling  surface  (Figure  35)  becomes  non-negligible  at  about  the  same  geometry  and 
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continues  to  grow  well  into  the  repulsion  wall,  indicating  that  coupling  may  continueto 


manifest  itself  throughout  the  potential  well.  While  the  DCT  and  coupling  surface 
provide  some  insight,  thetruevalue  ofthis  sample  calculation  won't  be  known  until 
these  surfaces  are  used  for  wave  packet  propagation. 

Werner's  Coupling  Angle 

While  we  have  used  Werner's  method  of  approximating  the  coupling  angle  [12] 
to  validate  our  COLUMBUS  calculations,  our  conclusions  also  speak  strongly  to  the  utility 
of  the  Werner  angle.  From  our  results,  we  may  conclude  that  the  Werner  method  has 
been  successfully  extended  to  apply  to  spin-orbit  MRCI  wavefunctions.  While  a  single 
point  calculation  of  K  He  required  on  the  order  of  24  CPU-hours  to  compute  the  DCT- 
derived  coupling  angle,  calculation  of  the  Werner  angle  required  approximately  half  of 
that  time.  This  work  has  shown  that,  depending  on  the  fidelity  one  requires,  it  may  be 
preferable  in  future  calculations  to  simply  compute  the  Werner  angle  rather  than  the 
DCTs  in  the  interest  of  resource  management.  This  recommendation  comes  with  the 
caveat  that  the  ease  of  calculation  of  the  Werner  angle  depends  on  the  complexity  of 
the  DRT  and  of  the  molecular  geometry.  For  sufficiently  complex  molecules,  this  may 
not  be  a  feasible  calculation,  and  the  interested  party  may  need  to  return  to  the  DCT- 
based  calculation. 
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Recommendations  for  Further  Work 


This  dissertation  is  not  an  end  in  and  of  itself,  but  has  provided  a  tool  for  many 
future  works. 

Spin-Orbit  Density  Matrix 

In  the  formalism  presented  here,  we  have  introduced  a  new  mathematical  tool, 
the  spin-orbit  density  matrix  (see,  e.g.,  equation  (185)).  While  the  standard  one- 
electron  density  matrix  exhibits  a  number  of  interesting  properties  (e.g.,  the  trace  is 
equal  to  the  total  number  of  electrons,  and  each  diagonal  element  gives  the 
approximate  electron  occupation  of  a  given  orbital),  we  have  not,  in  this  work,  explored 
the  full  meaning  and  properties  of  the  spin-orbit  density  matrix.  Such  an  exploration 
could  bring  new  understanding  and  utility  to  the  spin-orbit  density  matrix. 

Relativistic  Geometry  Optimization 

With  the  energy  gradient  method  in  place,  it  will  certainly  prove  enlightening  to 
compare  current  non-relativistic  optimized  geometries  to  those  calculated  with  the 
spin-orbit  Hamiltonian.  This  analysis  could  reveal  how  important  the  spin-orbit  effects 
are  to  the  geometries  of  simple  diatoms  to  complex  proteins. 

Wave  Packet  Propagation 

The  diabatic  surfaces  computed  herein  may  now  be  used  in  wave  packet 
propagation  methods  to  further  explore  the  dynamics  of  K  He.  The  diabatic  surfaces  of 
similar  systems,  such  as  Rb  He  and  Cs  He  can  also  be  calculated  with  similar  resources 
and  likewise  employed. 
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Analysis  of  Angular  DCTs 


Because  of  our  interest  in  the  radial  coupling  mechanism,  this  dissertation  has 
primarily  focused  on  the  radial  DCT,  which  is  calculated  as  the  z-component  of  the 
position  of  either  the  potassium  or  helium  atom  changes.  However,  the  DCT  is  truly  a 
multi-atom  gradient  with  a  different  component  for  each  of  the  six  coordinates  (in  the 
case  of  K  He).  Figure  38  shows  the  remaining  off-axis  DCTs  of  K  He.  The  first  property  to 
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Figure  38.  X-  and  Y-DCTs  of  K  He 


note  is  that  movement  of  the  helium  or  potassium  produces  different  couplings,  which 
is  to  be  expected  as  these  are  not  symmetric  atoms.  One  should  also  note  that,  despite 
the  supposed  symmetry  around  the  axis,  theX-  and  Y-  DCTs  for  helium  are  not  equal  as 
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the  geometry  approaches  the  repulsion  wall.  This  unexpected  phenom  enon  may  be  an 
artifact  of  forcing  the  use  of  C,  symmetry.  When  using  this  symmetry,  all  states  in  the 
SA-MCSCF  step  are  in  the  same  irrep,  and  must  be  unequally  weighted  in  order  to 
account  for  the  effects  of  orbital  resolution  on  the  density  matrices  (consider  equations 
(202)  and  (203)).  Whetherthis  compromise,  or  one  like  it,  is  the  cause  forthe 
difference  in  the  X-  and  Y-DCTs,  cannot  be  known  without  further  exploration;  however, 
given  the  agreement  with  the  Werner  coupling  angle,  the  answer  seems  to  have  no 
impact  on  the  fidelity  of  the  radial  DCT. 

The  Cartesian  DCTs  can  be  transformed  to  three  internal  DCTs  (radial  and  two 
angular)  and  three  center-of-mass  DCTs  by  means  of  a  Jacobi  transformation  [54],  The 

angular  DCTs,  shown  in  Figure  39,  are  of  the  form  ^IT^  and  are  in  fact  an  off- 

diagonal  element  of  the  nuclear  angular  momentum  operator,  J" .  This  equivalence  can 

provide  another  source  of  validation  of  the  DCTs,  if  J"  can  be  calculated  in  an 
alternative  manner.  Yarkony  [55]  contends  that  if  an  operator,  such  as  total  angular 
momentum,  commutes  with  the  Flamiltonian, 

W,i]  =  0  (272) 

then  it  is  truethat 

=  £,('P,|7|'P,)-<T,|7|T,)E,  (273) 
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Figure  39.  Angular  DCTs  of  K  He 


If  we  consider  that  this  total  operator  consists  of  both  nuclear  and  electronic 
components, 

J  =  Jn+Je  (274) 

then  equation  (273)  implies  that 

(xP/|7'!|xPy)  =  -(xP/|7e|xP/)  (275) 

That  is,  the  representation  of  the  nuclear  angular  momentum  operator  (or  angular 
DCTs)  is  equal  and  opposite  to  the  representation  of  the  electronic  angular  momentum 
operator.  In  non-relativistic  MRCI  calculations,  since  spin  is  not  considered,  the  total 
angular  momentum  is  L,  the  orbital  angular  momentum.  The  comparison  of  angular 
DCTs  to  Le  is  readily  made,  as  its  matrix  elements  are  quickly  calculated  by  COLUMBUS 
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from  the  available  orbital  angular  momentum  integrals.  However,  for  the  spin-orbit 


Hamiltonian,  in  which  orbital  angular  momentum  is  not  a  good  quantum  number, 

/V  /\  /V 

equation  (272)  does  not  hold  for  L,  but  rather  for  total  angular  momentum,  L  +  S  . 

Spin  integrals,  unlike  orbital  angular  momentum  integrals,  are  not  readily  available  in 
NWCHEM  or  COLUMBUS,  and  so  a  comparison  of  the  angular  DCTs  to  electronicangular 
momentum  is  not  possible  at  this  time.  If  NWCHEM's  integral  program  can  be  modified 
to  produce  the  spin  integrals,  and  COLUMBUS  to  produce  the  elements  of  the  spin 
operator,  then  another  source  of  validation  of  the  DCTs  would  be  available. 
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Appendix  A.  The  Hamiltonian 


This  appendix  provides  a  brief  introduction  to  the  energy  operator,  or 

Hamiltonian.  Energy  is  important  because  it  drives  the  dynamics  of  a  system,  the 

understanding  of  which  is  the  ultimate  goal  of  this  research. 

The  reason  for  the  preeminent  position  of  energy  eigenfunctions  in  the  theory  is, 
of  course,  that  they  form  the  stationary  states  of  isolated  systems  in  which 
energy  must  be  conserved. . . .  Moreover,  because  of  their  simple  time 
dependence. . .  they  can  be  used  to  build  up  a  description  of  the  time  evolution 
of  non-stationary  systems  as  well.  [56] 


Origin 

Formally,  the  Hamiltonian  is  derived  from  the  Lagrangian  via  the  Legendre 
transformation  [25]: 

H  ( q,p,t )  =  'Yjqipi-L(q,q,t )  (276) 

i 

where  q,  q,  and  p  are  generalized  position,  velocity,  and  momentum,  respectively,  and 
the  Lagrangian  is  defined  as  the  difference  of  kinetic  and  potential  energies: 

L  =  T-V  (277) 

However,  when  the  Lagrangian  of  a  system  can  be  expressed  as  a  sum  of  functions 
which  do  not  mix  degrees  of  generalized  velocities,  that  is, 

L  =  4  (q,  0  +  T,  (q,  q,t)  +  L2  (q,  q,t)  (278) 

where  L,  is  of  first  degree  in  q  and  L2  is  of  second  degree  in  q,  then  the  Hamiltonian 
is  just  the  total  energy: 
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H  —T  +  V 


(279) 


In  our  quantum  mechanical  treatment  of  nuclei  and  electrons,  the  kinetic  energy  will 
only  depend  upon  the  square  of  the  velocity  and  the  potential  energy  will  have  no 
velocity  dependence,  and  so  it  follows  that  the  Lagrangian  does  indeed  take  the  form  of 
equation  (278),  and  thus  the  Hamiltonian  will  take  the  form  of  equation  (279).  Although 
the  above  equations  are  classical  in  nature,  often  the  procedure  for  creating  quantum 
mechanical  formulae  is  to  merely  replace  the  classical  functions  with  operators  [26], 

This  will  be  the  case  forthe  Hamiltonian  operator,  which  will  takethe  form 

H=T  +  V  (280) 


Kinetic  and  Potential  Energy  Operators 

From  the  classical,  non-relativistic  analog,  the  kinetic  energy  operator  will  take 
on  the  form 


T  = 


2m 


(281) 


where  p  is  the  (linear)  momentum  operator.  To  extract  the  form  of  the  momentum 
operator  in  position-space,  first  consider  that  diffraction  experiments  indicate  that  a 
free  particle  with  a  known  momentum  has  a  wavefunction  of  the  form  [40] 


y/(r,r)  =  exp 


-( p-r-et ) 


(282) 


where  r  is  the  position,  s  is  the  energy,  and  ft  is  the  reduced  Planck's  constant.  It  is  a 
postulate  of  quantum  mechanics  that  the  operator  of  an  observable  quantity  is 
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hermitian,  and  thus  obeys  an  eigenvalue  equation  with  real  eigenvalues  [26];  in  the  case 


of  momentum  and  energy, 


py/  =  py/ 

Hi//  -  sy/ 

To  preserve  these  relationships,  the  momentum  operator  in  the  position-space 
representation  (as  opposed  to  the  momentum-  or  k-space  representation)  and  the 
energy  operator  take  the  form 


(283) 


p  — >  -ifN 

H^i  — 
dt 


(284) 


and  so  the  non-relativistic  Hamiltonian  will  always  contain  a  second  spatial  derivative. 

The  form  of  the  potential  energy  operator  will  depend  on  the  nature  of  the 
particular  system;  however,  we  can  say  generally  that  it  will  be  a  function  of  position 
but  not  of  velocity, 

V->V(r)  (285) 

(in  the  case  of  a  free  particle,  this  potential  is  zero).  Thus  we  can  cast  the  Hamiltonian 
as  a  solvable  differential  equation  in  position  space. 


H  =  -  —  V2  +V(r) 
2  m  V  ’ 


(286) 


The  Schrodinger  Equation 

The  Schrodinger  wave  equation  defines  the  relationship  between  the  time 
derivative  of  a  wavefunction  and  its  spatial  derivative.  For  a  single  particle,  it  is 
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(287) 


Hi// 


V  I 

-^-+y(F)  y/=iiy/ 

2m 


and  is  one  of  the  postulates  of  quantum  mechanics  [26],  In  equation  (284)  we  found 
that  the  energy  operator  is  expressed  as  ift,  which  also  appears  on  the  right-hand  side 
of  equation  (287)  with  eigenvalues  of  energy,  s  (see  equation  (283)).  Consequently, 
the  time-derivative  operator  can  be  replaced  by  its  eigenvalues,  and  the  Schrodinger 
equation  itself  takes  the  form  of  an  eigenvalue  equation, 


Hi y 


+  V(r)  y/ =  sy 


(288) 


This  is  the  time-independent  Schrodinger  equation  [26],  The  Hamiltonian  may  change 
from  system  to  system,  but  the  basic  eigenvalue  relationship  with  the  wavefunctions 
and  energy  will  remain.  It  is  this  eigenvalue  relationship  that  drives  the  mathematics  of 
quantum  chemistry. 
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Appendix  B.  The  Breit-Pauli  Hamiltonian 


The  Dirac  Equation 

Quantum  chemistry  relies  on  approximate  solutions  of  the  Schrodinger  equation, 
which  itself  can  be  considered  an  approximation.  If  wielded  carefully,  this  equation  will 
yield  most  of  the  information  about  the  energy  of  a  system;  however,  in  some  cases  it  is 
insufficient  for  the  accuracy  required.  By  assuming  the  kinetic  energy  takes  the  form  in 
equation  (281),  which  leads  to  the  Hamiltonian  in  equation  (8),  we  have  neglected  some 
of  the  less  significant,  relativistic  effects  one  finds  when  electron  speed  approaches  a 
non-negligible  fraction  of  the  speed  of  light. 

Close  scrutiny  of  the  Schrodinger  equation  ((288))  reveals  that  it  is  not  Lorentz 
invariant;  that  is,  it  does  not  treat  space  and  time  on  equal  footing.  A  more  suitable 
equation  for  exploring  relativistic  effects  is  the  Dirac  equation  fora  single  electron  [57] 


p  +—A  +  J3c 
c 


(289) 


where  A  is  the  vector  potential,  rp  is  the  scalar  electrostatic  potential,  a  is  a  three- 
component  vector  of  4x4  matrices,  and /?  is  a  single  4x4  matrix  [57],  Although  at  first 
glance  this  equation  appears  to  be  the  Schrodinger  equation  with  a  different 
Hamiltonian,  the  matrix  nature  of  the  equation  means  that  the  solutions,  \p ,  are  four- 
component  spinors  ratherthan  scalar  wavefunctions.  The  electrostatic  potential  may 
include  the  nuclear  potential,  and  the  vector  potential  has  been  included  in  the 
canonical  or  generalized  momentum  [58]  [59],  Since  the  momentum  only  appears  once, 
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the  Dirac  equation  has  only  a  first  spatial  derivative,  symmetric  to  the  first  time 
derivative.  Solutions  to  the  Dirac  equation  simultaneously  satisfy  the  simple  relativistic 
wave  equation,  the  Klein-Gordon  equation,  the  continuity  equation,  and  the  principle  of 
linear  superposition  [57], 


The  Breit  Equation 

When  introducing  additional  electrons,  an  electron  interaction  term  must  be 
added  to  equation  (289).  The  non-relativistic  two-electron  operator. 


1 

2 


(290) 


where  r  =  r  —r.  ,  assumes  instantaneous  Coulombic  forces  and  thus  is  not  Lorentz 

invariant.  A  fully  covariant  description  of  the  electron  interaction  lies  in  the  complicated 
Bethe-Salpeter  equation  [60],  A  truncation  of  the  expansion  of  that  equation  leads  to 
the  Breit  operator 


a,i  'ai  lJ  'ij  ij 


a.  -  a .+ 

1  j 


ro 


(291) 


which,  while  not  completely  Lorentz-invariant,  is  correct  to  the  order  of  a4 ,  where  a  is 
the  fine  structure  constant  (  =  cl)  [60],  Sincethis  operator  involves  the  direct  product 
of  the  a  matrices,  it  will  involve  16x16  matrices.  The  resulting  electronic  Hamiltonian 


(cf.  equation  (8))  is 
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a,*  '  ai 


The  potential  term  involving  the  nuclei  retains  its  non-relativistic  form  and  is  included  in 
the  electrostatic  potential  term  in  equation  (289): 


(293) 


r  . 

a,i  ' ai 


By  means  of  a  Foldy-Wouthuysen  transform,  the  Hamiltonian  is  reduced  into  this  final 
form,  the  Breit-Pauli  Hamiltonian  [60]: 


"  2  ft,  ftj  1 

'“?c  ”w+T~i7+l1 


(<T  ’  B, ) ~  <x,  •  (^  x  %  -  E,  x  ki  ) 


8c 


-  2(' 


(294) 


1 


4c2r3 
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J  2c2/;.3  L 

v 


where  ;r(.  is  the  gauge-covariant  momentum  of  electron  /  [32],  is  the  total 


electric(magnetic)  field  at  electron  /,  and  cr.  is  the  spin  of  the  ith  electron.  Tinkham 
justifies  this  simplification  to  the  Breit-Pauli  Hamiltonian  by  the  relatively  low  energies 
addressed  in  quantum  chemistry  [56], 

This  Hamiltonian  can  be  sorted  into  seven  operators: 

1.  The  non-relativistic  Hamiltonian 


H, 


ru 


(295) 
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2.  The  mass-velocity  correction  term 


».=-! 


Pi_ 

8c2 


(296) 


3.  The  orbit-orbit  term 


u,= 


Z 

*.7 


2c- 


\Pi-pj)+ — 


)0v*0 


4.  The  spin-orbit  term 


(297) 


H3=~K 

3  2c2 


X^'t^xA)+X2^ 


IJ  — 

V  y  7 


5.  The  Darwin  term 


»4=Z^M,) 


(298) 


(299) 


6.  The  spin-spin  term 


2;r 1  _ 
3  ^ 
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(cr  -T.)( 
r 

'j 


°i 


- 


(300) 


7.  The  external  magnetic  interaction  term 

Ht=i'Ll(vrB,)  +  (A,-p,)  (301) 

i 

The  electron  rest-mass  energy,  c2,  also  appears  in  the  Hamiltonian,  but  will  be  ignored 
since  it  serves  only  as  an  energy  offset.  Although  the  Dirac  equation  acts  on  four- 
component  spinors,  the  above  seven  terms  have  scalar  functions  as  their  eigenvectors. 
"[T]here  is  really  no  longer  any  serious  difficulty  in  principle  in  writing  down  a  basic 
Hamiltonian  which  is  accurate  enough  for  all  practical  purposes  in  extranuclear  physics" 
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[56],  Thus,  they  can  be  appended  to  the  non-relativistic  Hamiltonian  in  the  Schrodinger 


equation,  leading  to  more  accurate  eigenvectors  and  eigenvalues. 

Fine  Structure  and  the  Spin-Orbit  Operator 

The  orders  of  magnitude  of  the  Hx,  H3,  and  H4  terms  are  on  the  order  of  a 2 
times  smallerthan  the  non-relativistic  Hamiltonian;  these  are  grouped  together  as  the 
fine  structure  terms  [61].  The  others  are  the  hyperfine  structure  terms ,  and  are  at  least 
an  order  of  magnitude  smaller  than  the  fine  structure  terms  [61].  Forthis  reason,  one 
can  disregard  the  hyperfine  structure  terms  and  still  get  a  much-improved  solution 
compared  to  the  non-relativistic  Hamiltonian. 

Let  us  focus  on  the  spin-orbit  term,  H3 ,  as  the  other  fine  structure  terms,  Hx 
and  H4,  have  little  effect  in  the  valence  region.  By  dissecting  the  electric  field  term  in 
equation  (298),  we  can  see  the  spin-orbit  term  separates  into  one-  and  two-  electron 
pieces.  The  field  is  due  both  from  the  nuclei  and  the  electrons: 


r.  + 


I 


—  r. 
r3  11 

V 


(302) 


(where  no  external  field  is  assumed)  and  so  the  cross  product  Et  x  p.  is 


Ei*Pi  =-Z%(4xa)  +  Z4(^xa)  (303) 

a  ria  j  rij 

We  recognize  (ria  x p,)as  the  orbital  angular  momentum  of  the  ith  electron  with  respect 
to  the  a'h  nucleus,  and  (r.  x  p.)  as  the  angular  momentum  of  the  i'h  electron  with 
respect  to  the  jth  electron: 
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(304) 


(4  xA)  =  4 

Making  these  substitutions,  equation  (298)  takes  the  form 


#,=■ 


2c2 


-V  Z 


— fcr.  L,„  +  >  —  cr. 


'■V  ij 


(E  -2L-) 

V  ij  }' ) 


(305) 


Equation  (305),  more  so  than  equation  (298),  shows  the  separate  one-  and  two-electron 
spin-orbit  terms.  We  can  notate  these  as 


a  'ia 

4°(u)=^E^-(4-24) 

z(  y  'ij 


(306) 


The  one-electron  term  grows  rapidly  with  respect  to  nuclear  size  as  compared  to  the 
two-electron  term,  and  so  the  latter  is  less  important  for  larger  elements  [62],  In  any 
case,  the  effects  of  the  two-electron  term  may  be  effectively  included  in  the  AREP.  The 
revised  electronic  Hamiltonian  is 


H  = 


i  Z  Uj 


ri~rj\ 


-z 


r  -r 


a,  ■  L. 


(again,  cf  equation  (8)). 


(307) 
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Appendix  C.  Creation  and  Annihilation  Operators 

In  this  appendix,  we  discuss  the  creation  and  annihilation  operators,  used  to 
second-quantize  the  Hamiltonian,  in  greater  depth.  While  not  a  proof,  this  section 
employs  a  simple  example  involving  small  Slater  determinants  that  the  reader  should 
find  helpful  in  demonstrating  their  properties.  In  this  section,  we  will  use  kets  to 
represent  both  Slater  determinants  and  functions— to  differentiate,  we  leave  electron 
dependence  off  of  the  determinant  notation. 


Annihilation 

For  a  given  function,  we  recognize  that 

I  Xa  (0>®  I  Xp{i))®  I  (k))®...  =1  {i)Xp{j)Xy  (k)..) 
Considerthe  action  of  theoperator 

4nYj(X i(0' 

i 

on  an  example  Slater  determinant, 


\XxXiX3)  =  ^ 

W1)) 

ki(2)) 

^2(1)) 
^(2))  | 

1  ^ ^ 
cT 

_  cn 

_ki(3)> 

^2(3)) 

X3 (3))J 

This  is  equal  to 


(308) 


(309) 


(310) 


(!) 1  +(Xi  (2)  I  +<Xi  (3)  l)  (I  Xl  (1)  X2  (2)x3  (3))- 1  Xi  (1)^3  (2)  J2  (3))  + 

I  X3  (!)  Xx  (2)  J2  (3))- 1 (1) Xx  (2)  X3  (3))+ 1 J2  (1)  (2)Xi  (3))- 1 11)  ^2  (2) Xi  (3)» 
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A  (%\  (i)  I  bra  will  only  create  a  non-zero  inner  product  with  a  ket  that  includes  I  (i))  in 

its  direct  product.  The  resultant  inner  product  will  integrate  out  the  i'h  electronic 
variable  since,  by  orthonormality,  we  assume 

fc(0l^(0)=<^  (312) 

The  resultant  sum  of  functions  is 

iO  ^2  (2 )z3  (3)>- 1  X3  (2)z2  (3)>+ 1  (1  )%2  (3)>)-  3 

^(1  X 2  (1)^3  (3)>+  I  (1)^3  (2))-  I  ^3  (1)^2  (2)>) 

While  this  does  not  appear  to  be  a  Slater  determinant  at  first  glance,  when  each  pair  of 
electrons  is  considered  alone,  we  see  that  there  are  in  fact  three  determinants: 

I  X2X3)+ 1  X3X2)+ 1  X2X3)  (314) 

Since  I X2X3)  =  ~\  X3X2)  >  equation  (314)  reduces  to  the  determinant  I  X2X3) .  Thus  the 
operator  in  (309)  has  annihilated  the  X\  spin-orbital  from  the  determinant.  For  that 
reason,  we  can  equate  it  to  an  annihilation  operator  [40]: 

^N^ix i(0|s«i  (315) 

i 

with  the  action 

a,\  X1X2X3)  X2X3)  (316) 

If  we  were  to  split  up  the  spin-orbital  into  its  spatial  and  spin  components,  we  define 

(31?) 
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Creation 


Next  consider  the  action  of  the  operator 

(318> 

on  the  Slater  determinant  I  X2X3)  ■  The  operator,  by  adding  an  orbital  to  the  product, 
must  effectively  also  add  a  third  electron.  Thus,  in  three-electron  space,  there  are  three 
forms  of  the  determinant  \  X2X3)  that  must  be  considered  in  the  product: 

i 0*2  (1)X3  (2))-  lz3  (1)Z2  (2)>) 

-i('^(l)Z3(3)HZ3(l)^(3)))  (319) 

^(l  ^2  (2)^3  (3)>- 1  (2)^2  (3))) 

(Note  that  the  second  determinant  was  multiplied  by  -1.  This  is  because  of  the  anti¬ 
symmetric  nature  of  electrons,  which  requires  that  the  wavefunction,  upon  exchange  of 
electrons,  must  have  opposite  sign.  The  sign  change  is  not  present  in  the  third 
determinant  because  two  electrons  were  exchanged,  thus  cancelling  the  sign  change 
[26].)  The  direct  product  with  the  operator  is  then 

i(U,(l)>+U,(2)>+U,(3)>)®[(lJl(l)j3(2)>-U,(l)3rJ(2)>)- 
0^(1)^(3)>-I^(1)^(3)>)  +  (I^(2)^(3)>-I^(2)^(3)>)] 

A  I  X\  (*))  ket  will  only  create  a  non-zero  direct  product  with  a  ket  that  does  not  already 

have  a  function  of  the  f‘  electronic  variable.  The  resultant  function  is 

^[(!  ^2  (1)*3  (2) Xx  (3))- 1  X3  (1)^2  (2) Zi  (3))) 

-(l  J2  (l)Xi  (2)X3  (3))- 1  X3  (l)^i  (2)Xi  (3)))  (321) 

+  (l  Zi  (0  X 2  (2) X3  (3)>- 1  Xi  (!)  (2)  2f2  (3)))] 
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which  is  the  determinant  I  Z1Z2Z3)  ■  Thus  the  operator  in  (318)  has  created  the  Z\  spin- 
orbital  in  the  determinant.  For  that  reason,  we  can  equate  it  to  a  creation  operator  [40]: 

,322) 

with  the  action 

I  Z2X2)  =l  Z1Z2Z3)  (323) 

As  in  the  case  for  the  annihilation  operator,  we  can  split  the  spin-orbital  into  spatial  and 
spin  components: 

(324) 

The  creation  and  annihilation  operators  obey  the  anticommutation  relationship  [40]: 

{a*  ,d  >  =  a*  a  +  d  d t  =1.  (325) 

These  are  the  creation  and  annihilation  operators  introduced  in  the  main  body  of  the 
text. 
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Appendix  D.  The  Spatial  Orbital  Basis 


In  order  to  represent  the  Hamiltonian,  the  first  task  is  to  build  a  suitable  basis  for 
a  single  electron;  from  that  set,  a  basis  for  the  multi-electron  wavefunctions  can  be  built 
as  a  direct  product  space.  This  one-electron  basis  must  meet  three  criteria: 

1.  It  must  be  orthonormal 

2.  It  must  be  easy  to  integrate  (for  computational  savings) 

3.  It  must  predict  the  lowest  possible  energy  (as  a  consequence  of  the  variation 
principle) 

Orthonormality 

For  asingle  atom,  thedirect  product  ofspherical  harmonics  in  0and  rp ,  with 
associated  Laguerre  polynomials  in  r  provides  a  set  of  orthonormal  orbitals  with  which 
to  construct  wave  functions  [28],  These  orbitals  are,  in  fact,  the  eigenfunctions  of  the 
hydrogen  atom  Hamiltonian  [26]  and  are  known  as  Slater  Orbitals  [20];  however,  for 
polyatomic  molecules,  hydrogenic  basis  functions  located  at  different  atomic  centers 
are  not  orthogonal,  having  an  overlap  matrix  defined  by 

s^=(^k.')  <326) 

where  the  %k  are  the  orbitals.  Such  a  basis,  which  we  will  refer  to  as  an  atom-centered 
basis  or  atomic  orbitals,  must  be  transformed  non-unitarily  into  a  non-localized  set  of 
orthonormal  Molecular  Orbitals  (MO),  often  referred  to  as  a  Linear  Combination  of 
Atomic  Orbitals  (LCAO)  [20],  These  MOs  are  defined  as 
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(327) 


h)=HcMiM 


Since  the  molecular  orbitals  are  required  to  be  orthonormal,  a  transformation  will 
reduce  the  overlap  matrix  to  the  identity  matrix. 

Integratability 

Slater  orbitals,  while  more  physically  accurate,  are  unwieldy  for  computation. 
For  example,  the  IS  Slater  orbital  is 


(328) 


\n  J 


which  exhibits  cusp  behavior  at  the  origin.  Figure  40  shows  the  square  root  of  the 
probability  of  finding  an  electron  as  a  function  of  R ,  in  atomic  units.  Integration  against 
other  Slater  orbitals  at  other  origins  proves  to  be  computationally  intensive.  For  this 
reason,  an  alternative  set  of  Gaussian  orbitals  is  often  chosen  instead.  The  pivotal 
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Figure  40.  IS  Slater  Orbital,  £=1.0 
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difference  between  the  Slater  orbital  and  the  Gaussian  orbital,  which  has  the  form 


(329) 


V  n  7  v  ' 

is  that  in  the  latter  the  argument  is  squared  in  the  exponent.  When  two  Gaussian 
functions  are  integrated,  the  result  is  a  Gaussian,  regardless  of  their  origins,  making  the 
task  of  integration  much  easier.  Figure  41  compares  a  Gaussian  orbital  to  the  previous 
Slater  orbital.  The  Gaussian  is  thicker  towards  the  origin  and  falls  off  faster,  while  it 
exhibits  no  cusp.  Although  it  is  not  a  very  accurate  representation  of  the  Slater  orbital, 
the  ability  to  develop  efficient  integration  algorithms  with  Gaussian  orbitals  means  that 
more  of  these  can  be  used.  Figure  42  shows  a  linear  combination  of  three  Gaussian 
orbitals  against  the  original  Slater  orbital.  The  figure  shows  that,  although  not  quite  a 
match  at  the  cusp,  this  orbital  is  much  closer  to  the  Slater  orbital,  and  is  still  more  easily 
integrated.  A  basis  created  from  such  functions  is  labeled  a  STO-3G  basis.  STO  indicates 
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Figure  41.  Gaussian  orbital,  a=0.6 
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that  it  is  a  Slater-Type  Orbital  (as  opposed  to  a  Slater  orbital)  and  3G  indicates  the 
orbitals  are  each  made  from  three  Gaussians  [20],  Although  the  HF  method  itself  can  be 
applied  to  either  Slater  orbitals  or  STOs,  the  integration  advantage  of  the  STOs  makes 
them  a  staple  in  computational  chemistry  packages. 

Lowest  Possible  Energy 

The  variation  principle  [20]  states  that,  for  an  approximate  eigenvector^) ,  the 
eigenvalue  will  be  related  to  the  true  energy  by 


{y/\H\y/)>s0 


(330) 


Thus,  given  several  different  bases,  the  one  that  produces  the  lowest  eigenvalue  is  the 
best  choice. 
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Figure  42.  Linear  combination  of  three  Gaussian  orbitals  [20] 
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Appendix  E.  Self  Consistent  Field  Method 


The  Fock  Operator 

For  our  purposes,  the  Hartree-Fock  method  results  in  the  best  possible  ground 
state  of  a  molecule  given  a  molecular  orbital  basis  derived  from  a  single  Slater 
determinant.  While  we  will  ultimately  be  interested  in  optimizing  the  basis  for  excited 
states,  this  method  nevertheless  provides  a  first  guess  from  which  a  more  suitable  basis 
may  be  formed.  From  the  preceding  section,  we  have  found  that  the  orbitals  must 
satisfy  orthonormality, 

=  (33D 

and  are  defined  by  the  variation  principle.  Using  a  single  Slater  determinant  for  the 
ansatz  y/  leads  to  the  constraint  equation  [20] 


jdjfi(j)  — 


^(0-ZJ 

j*i 


r.J 


4ifi(j)ta(j)—  s K  (0  =  (0  <332) 


where  h  is  the  single-electron  operator  from  equation  (13).  The  integral  in  the  second 
term  in  equation  (332)  is  called  the  Coulomb  operator,  Jb  (j)  and  represents  the 

average  instantaneous  coulomb  potential  between  electron  i  and  electron  j  .  The  third 
term  in  equation  (332)  looks  almost  exactly  the  same,  except  that  the  electrons  have 
been  exchanged  between  the  two  orbitals.  For  this  reason,  this  integral  is  known  as  the 
exchange  operator,  and  has  the  following  property: 


Kb(i)</>a{i)  = 


(333) 
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The  exchange  operator  is  a  result  of  the  antisymmetric  nature  of  electrons  (due  to  the 


Pauli  exclusion  principle),  and  doesn't  have  a  simple  classical  analog.  The  form  of  the 
operator  in  equation  (333)  allows  us  to  rewrite  equation  (332)  as 


‘(O+EMO-EMO 


J*l 


J** 


t  (0 = £A  (0 


(334) 


The  operator  on  the  left  is  called  the  Fock  operator,  f  (/)  [20],  The  orbitals  must  obey 
the  pseudo-eigenvalue  equation 

f\ti)  =  £iW  (335) 

which  are  the  Hartree-Fock  equations.  This  is  not  a  true  eigenvalue  equation  because 
the  operator  /  depends  upon  the  eigenvectors  |$),and  it  must  be  solved  iteratively. 

That  is,  a  guess  for  the  eigenvectors  is  proffered,  from  which  the  Fock  operator  may  be 
calculated,  from  which  new  eigenvectors  can  be  calculated,  and  so  on  until  the 
eigenvectors  reach  self-consistency. 


Pople-Nebset  Equations 

The  HF  method  is  further  subdivided  into  a  number  of  different  cases.  We  will 
not  review  them  all,  but  instead  will  present  one  of  the  more  general  of  these,the 
unrestricted  open-shell  Hartree-Fock  method  (UH F),  as  an  example. 

Unrestricted  refers  to  the  fact  that  spin-up  electrons  will  be  allowed  a  different 
set  of  spatial  orbitals  than  the  spin-down  electrons.  Thus  equation  (327)  is  split  into  two 
separate  equations  for  spin-up  and  spin-down: 
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There  are  also  two  Fock  operators, 


/“  =(a\f\a) 

ff  =  {P\f\P) 

Plugging  these  definitions  into  the  HF  equations  leads  to,  for  example. 


(337) 


(338) 


which,  when  integrated  against  1%V  | ,  yields 


(339) 


where  S  is  the  overlap  matrix  (see  equation  (326))  and  F  is  the  Fock  matrix  whose 
elements  are  defined  as 


K=(xXf\xM)  (340) 

with  symmetric  equations  for  spin-down  (note  this  is  not  the  same  second-quantized 
Fock  matrix  used  in  the  formalism  in  the  main  body  of  the  text).  Equation  (339)  is 
equivalent  to  the  matrix  equation 

ForCor=SarCVr  (341) 

which  are  known  as  the  Pople-Nebset  equations  [20],  By  transforming  this  equation  into 
a  basis  in  which  the  overlap  matrix  becomes  the  identity  matrix,  the  Pople-Nebset 


equations  become  pseudo-eigenvalue  equations.  This  transformation  is  effected  as 
follows: 

1.  Since  S  is  hermitian,  diagonalize  it  by  the  transform 

s  =  UfSU  (342) 

2.  Create  the  matrix  s  a  diagonal  matrix  whose  entries  arethe  inverse  of  the 
positive  square  roots  of  the  eigenvalues  of  s . 

3.  Define  the  transformation 

X  =  Us-^Ut  (343) 

4.  Using  this  transformation,  create  the  matrices 


C" '  =  X_1C“ 

F" '  =  XfF"X 

along  with  their  spin-down  counterparts,  to  create  the  equations 


(344) 


F«  ' c“ '  =  Ca  ' «?“  (345) 

Both  F3  and  F/J  depend  on  C"and  C/y,so  the  two  equations  are  not  independent. 
Thus,  in  order  to  solve  these  equations,  one  must: 

1.  Make  an  initial  guess  at  F“  and  F^ 

2.  Transform  F“ 'and  F^' 

3.  Calculate  their  eigenvectors  and  eigenvalues;  equate  these  to  C“ Cp\  ^“and 

£P 

4.  Transform  these  back  to  C"  and  Cp ,  use  these  to  calculate  a  new  F“  and  F/J 
and  so  forth,  until  self-consistency  is  achieved. 
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When  the  single-electron  orbitals  have  been  determined,  the  final  step  of  the 


SCF  method  is  to  create  an  antisymmetric  Slater  determinant  from  the  lowest  N 
orbitals: 

h  )=k&--4)  (346) 

This  is  the  electronic  ground  state,  used  to  determine  the  ground-state  energy 

= W^ko)  (34?) 

As  a  final  note,  since  C  is  a  square  matrix,  there  will  be  as  many  MOs  created  as 
AOs  used,  regardless  of  the  number  of  electrons  involved.  Because  the  MOs  not 
occupied  by  an  electron  don't  affect  the  energy,  those  unoccupied  orbitals  are 
orthonormalized  but  not  energy -optimized  by  this  process.  If  one  is  interested  only  in 
the  ground  state,  these  extra  orbitals  are  of  little  value  anyway;  however,  if  one  is 
interested  in  excited  states,  these  orbitals  are  not  ideal.  The  MCSCF  method,  to  be 
briefly  discussed  next,  addresses  this  concern  by  further  performing  the  SCF  method  to 
higher  states.  Although  the  SCF  MOs  can  be  used  directly  in  a  subsequent  Cl  calculation, 
an  intermediate  MCSCF  step  can  reduce  the  effort  of  the  Cl  by  identifying  and 
eliminating  those  orbitals  that  will  have  less  impact  on  the  Cl  result. 

Multi-Configuration  Self-Consistent  Field 

MCSCF  is  a  complicated  but  useful  generalization  of  the  SCF  procedure.  While 
the  SCF  method  optimizes  the  molecular  orbitals  by  adjusting  the  SCF  coefficients  (see 
equation  (327)),  and  the  Cl  method  optimizes  the  Cl  coefficients  of  the  determinants  to 
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create  the  eigenfunctions  of  the  Hamiltonian  (see  equation  (351)),  the  MCSCF  method 


combines  both  of  these  optimizations,  but  usually  on  a  smaller  scale  than  the  Cl: 

1.  The  ground-state  orbitals  are  optimized 

2.  Several  determinants  are  constructed  (not  just  the  ground  state) 

3.  A  linear  combination  of  determinants  is  optimized  forthe  ground-state  energy 

4.  New  orbitals  are  optimized  based  on  the  Cl  coefficients 

The  process  is  continued  until  self-consistency  of  the  involved  orbitals  is  achieved.  This 
procedure  specifically  optimizes  the  orbitals  for  the  ground  state;  the  method  can  also 
be  adapted  to  optimize  the  orbitals  for  some  excited  state.  Forthe  construction  of 
derivative  coupling  terms,  there  are  several  (but  commonly  two)  states  of  importance. 
Rather  than  optimizing  to  any  single  state,  the  orbitals  are  constructed  to  have  the  best 
average  optimization  for  all  states  involved.  This  method  is  called  the  State-Averaged 
MCSCF  (SA-MCSCF)  [3]. 

The  scope  of  this  project  does  not  involve  the  modification  of  the  MCSCF 
procedure,  and  only  relies  upon  it  inasmuch  as  it  prepares  orbitals  for  a  multi -electron 
basis;  hence,  we  will  not  give  more  than  this  rudimentary  explanation  here. 
Nevertheless,  the  interested  reader  may  wish  to  refer  to  Shepard's  paper,  which 
thoroughly  outlines  the  subject  [35], 
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Appendix  F.  Configuration  Interaction 

The  method  of  Configuration  Interaction  (Cl)  is  an  improvement  upon  the  SCF 
method  [20],  It  is  also  the  principle  computational  tool  with  which  we  calculate 
eigenfunctions  for  the  derivative  coupling  terms.  The  goal  of  Cl  is  to  solve  the  matrix 
equation 

Hc'=f'c7  (348) 

where  His  the  Hamiltonian  matrix,  s1  represents  an  energy  eigenvalue  corresponding 
to  the  eigenvectors  c1 ,  which  fulfill  the  relationship 

k/)  =  Zc/k)  (349) 

i 

While  the  SCF  method  used  the  lowest  N  orbitals  to  create  a  single  Slater  determinant, 
disregarding  higher  orbitals.  Cl  uses  some  or  all  of  the  orbitals  to  create  a  family  of 
Slater  determinants  in  which  to  represent  the  Hamiltonian.  The  ground  state  of  this 
Hamiltonian  will  no  longer  be  a  single  determinant  (see  equation  (346)),  but  a  linear 
combination  of  determinants.  This  Hamiltonian  yields  two  advantages  over  the  SCF 
Hamiltonian:  first,  the  ground  state  will  be  at  least  as  good  as  the  SCF  ground  state, 
leading  to  an  equal  or  lower  (and  hence  better)  ground-state  energy;  and  second,  the 
Hamiltonian  allows  for  excited  states. 

Given  a  system  of  n  spin-orbitals  and  N  electrons,  one  can  make  a  total  of 
Slater  determinants.  Each  determinant  may  be  constructed  by  replacing  an  occupied 
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ground-state  MCSCF  orbitals  with  one  of  the  unoccupied  (or  virtual)  MCSCF  orbitals.  For 
instance,  if  the  spin-orbital  |^,)  is  included  in  the  ground  state,  one  can  construct  all  the 
determinants  that  replace  \<j>^  with  an  orbital  not  included  in  the  ground  state,  \<j>  \ . 

This  is  a  singly-excited  determinant,  denoted  by  | (p'a  ^ .  There  are  N(n  —  N )  such  singly- 
excited  determinants  [20],  Similarly,  we  can  remove  two  spin-orbitals  at  a  time  and 
replace  them,  constructing  the  family  of  doubly-excited  determinants,  |  <prJb)>  of  which 

(N'\(n-N'\ 

there  are  .  Clearly  one  can  continue  constructing  determinants  in  this 

v2A  2 

fashion  until  all  combinations  are  exhausted.  Representing  in  such  an  N  -electron  basis 
gives  what  is  called  a  full  Cl  Hamiltonian: 

'  («,| «k)  (%\h\<p:)-  («,l«| <p1) 

v  •  •  • 

(If  n  =  N ,  the  Cl  space  reduces  to  a  single  determinant.)  Because  the  Flamiltonian  is 
hermitian,  only  the  upper  triangle  must  be  calculated.  Furthermore,  orthogonality 
restrictions  dictate  that  a  number  of  the  elements  will  be  zero  [20],  Nevertheless,  the 
size  of  this  matrix  grows  rapidily  with  the  number  of  basis  orbitals  and  electrons.  For 
this  reason,  larger  molecules  almost  never  receive  a  full  Cl  treatment.  Commonly,  the 
Flamiltonian  is  restricted  to  be  represented  only  in  the  ground  state,  singly-excited 
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determinants,  and  doubly-excited  determinants;  this  treatment  is  aptly  named  the 
Singles  and  Doubles  Cl  (SDCI  or  CISD)  [3],  This  method  calculates  even  lower  energies  if 
it  is  preceded  by  an  MCSCF  step. 

It  should  be  noted  that,  rather  than  using  bare  Slater  determinants,  the  CSFs  of 
the  Gel'fand-Tsetlin  basis  used  in  COLUMBUS  employ  a  spin-adapted  linear  combination 
of  Slater  determinants;  this  does  not  change  the  general  concept  discussed  here. 

Once  the  elements  of  the  Hamiltonian  are  calculated,  its  eigenvectors  are 
constructed  from  the  Slater  determinants, 

k')=4'ko)+Z(e)V.')+Z(s)>:>  (35i) 

ar  abrs 

by  a  matrix  diagonalization  or  similar  procedure.  These  are  the  wavefunctions  that  will 
be  used  to  calculate  derivative  coupling  terms. 

The  interested  reader  may  wish  to  reference  Shavitt's  more  complete  review  of 
the  Cl  method  [63], 

Multi-Reference  Cl 

COLUMBUS  uses  Multi-Reference  Singles  and  Doubles  Cl  (MRSDCI  or  MR-CISD). 
This  method  differs  from  the  standard  SDCI  in  that,  rather  than  taking  single-  and 
double-  excitations  only  from  the  ground  state  CSF,  these  excitations  are  taken  from  a 
set  of  CSFs,  known  as  the  reference  CSFs  determined  from  the  MCSCF  step.  In  cases 
where  information  is  needed  about  excited  states,  especially  nearly-degenerate  states, 
MR  methods  are  successful  where  single-reference  methods  may  fail  [9], 
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Direct  Cl 


The  Hamiltonian  matrices  calculated  in  computational  chemistry  can  have 
dimensions  on  the  order  of  hundreds  of  millions  of  elements  in  practical  calculations; 
however,  usually  it  is  only  the  lowest  handful  of  eigenvalues  and  eigenvectors  which  are 
needed.  Such  a  requirement  hardly  justifies  the  diagonalization  of  the  entire  matrix 
which,  for  an  nxn  Hamiltonian  can  require  on  the  order  of  n3operations  [64],  which  is 
too  costly  when  only  a  few  eigenvalues  are  required.  Furthermore,  calculations  in 
normal  diagonalization  routines  require  the  whole  Hamiltonian  be  used,  which  is  a 
fierce  burden  on  computer  memory.  Nesbet  proposed  a  method  [64]  which  Shavitt 
modified  [65],  now  known  as  Direct-CI,  which  uses  only  a  small  portion  of  the 
Hamiltonian  for  each  calculation,  producing  one  eigenvalue  and  eigenvector  at  a  time, 
so  that  the  effort  required  is  proportional  to  the  number  of  eigenvectors  sought.  This 
method  attempts  to  solve  for  the  eigenvalues  and  eigenvectors  of  equation  (348) 
iteratively,  without  solving  all  equations  at  once.  The  Davidson  method  [66],  itself  built 
upon  an  earlier  method  proposed  by  Roos  [67],  is  the  basis  for  the  current  direct  Cl 
method  used  in  COLUMBUS  [68], 
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Appendix  G.  The  Non-Relativistic  Energy  Gradient 


This  appendix  outlines  Shepard's  formalism  for  obtaining  the  non-relativistic 
energy  gradients  [3].  It  is  intended  to  be  a  quick  overview  for  readers  of  this  paper.  For 
a  more  in-depth  understanding,  please  see  Shepard's  referenced  paper. 

The  energy  gradient,  in  the  resolved  basis, 

s(Ry  =I^\R)\h[z\r)x\yY]{R))  (352) 

can  be  represented  in  the  [5]  basis  as: 

M Z]{R)  hW(r)xLW(r)\  = 

r  n,  (353) 

(y/f](i?)|  exp(-Z(^))exp(-^(/?))//[s](^)exp(^(i?))exp(z(^))  |^JS](/?)) 

When  the  exponentials  are  expanded  as  a  series,  they  form  commutators  with  the 
Hamiltonian  according  to  the  Baker-Cam pbell-Hausdorff  theorem: 

(<^z]  (tf)|tf[z]  (fl)A|^|z]  (R))  =  (ys\s]  (/?)|^[s](/?)Jr|^s](/?)) 

+  (vFp](/?)|[wIs](i?),Z(i?)J|^js](/?))  (354) 

+(^;s]  (/?)|[^[s]  |^;s]  (^))+- 

where  higher-order  terms  are  ignored  in  this  derivation  [3],  Since  each  of  the  rotation 
operators  is  unitary,  it  can  be  expressed  as 

exp  (A  (/?))  =  exp  £ara  (/?)(£„  ~Esr)  (355) 

(cf.  equation(115)).  Here,  we  have  combined  each  generator  with  its  transpose  so  that 
the  parameters  are  unrestricted.  (This  adjustment  was  primarily  made  to  coincide  with 
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Shepard's  notation  [3],  but  the  concept  is  unchanged.)  Thus  equation  (354)  can  be 
expressed  as 

(r!zH«)|HlzH«)’|r!zl('!))  =  (f'!s|(^)|H|sl('!)Vfl(^)) 

+  MS|(K)|[ZMR)[»|5l(tf)A-£jl  liT’W)  (356) 

V  rs  J 

+K^)|fZMfl)[»|sl(*)A, -£,,]]  WS1M) 

V  rs  J 

The  following  notational  substitutions  are  made  for  simplicity: 

X^W[h[s|WT„-£,.]=?W'/.S(«) 

(357) 

TJK,(R)[Hls](R).E„-E„]  =  k(R)-f“(R) 

r<s 

where  fCI  is  the  orbital  gradiant  vector,  and  the  dot  product  represents  a  sum  of 
operators  rather  than  a  scalar  quantity.  (There  is,  in  fact,  only  one  such  vector; 
however,  since  rotations  are  partitioned  into  redundant  and  essential  classes,  the 
orbital  gradient  vector  can  be  so  partitioned  as  well.)  Then  equation  (356)  is  simplified 
as 

(^Jz]  (/?)|^[z]  (/?)"  |^z]  (/?))  =  (/?)|^[s]  (/?)"  |^]  (i?)) 

(358) 

+{z(R)-m(R)]  +{k(R)-f^(R)] 
which,  expanded  and  evaluated  at  R0,  is  equal  to 

(r!zl  ( *0 )  I  ( *0 ) '  I  vT'  ( *0 ))  =  (kT1  ( R» )  I  ( *0  y '  I  (T1  ( K )} 

+(hkT ■ /:: w)  +{H^y-f:::w)  ossi 

+ ( z  ( )•/.*(  k  y )  +  (*  (  k  )■  iz  (  k  y ) 
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At  the  reference  geometry,  \y/a  {R))  =  \Vci  (^0))  anc*  so  ^(^0)  =  Z(^0)  =  0'  anc*  so 

their  parameters  must  be  zero.  This  leads  to  the  elimination  of  two  terms  in  the  above 
equation 


+(?(«„)'•  7r  K  ))  +  (k  (R,  )'•/.“  ( )) 


(360) 


Considerthe  first  term,  expanding  the  Hamiltonian  in  its  second-quantized  form: 


r!slK)|tfwK)' |ris| (R„j)  =  (rl!| (-Ro)|S(C (K)J  K 

rs 

+iE(sS2,  W)  er,n,|^!s]K) 

rstu 

=E(€1K))>^M£"K1W 

rs 

+il(*2w)‘(^w|«.>j,|w 


(361) 


We  use  the  density  matrices  (see  equation  (34))  to  simplify  equation  (361): 


r!s|  M\h  W  |r!s|  («„)) = Z(C  M)‘  oiI](«o) +iZ(«2.  W)'  42  M 


=  Tr 


(/1|s|(fi„))'.D|s|(^)l  +  W,{(*|s|(fi„))'.rf|s|(fi„) 


-,(362) 


where  the  nuclear  derivatives  now  only  act  upon  the  coefficient  matrices,  as  the 
generators,  and  hence  density  matrices,  are  geometry-independent. 

Now  consider  the  second  term  in  equation  (360).  The  orbital  gradient  vector  is 
non-zero  as  a  consequence  of  the  Cl  step,  which  optimizes  the  CSF  coefficients  but  not 
the  orbitals  [3],  For  an  MCSCF  step,  or  a  full  Cl,  the  orbital  gradient  vector  is  zero; 
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however,  forthe  sake  of  generality  we  will  assumethat  it  is  non-zero.  Its  form,  which  is 
derived  elsewhere  [5]  [35],  is 


u 


+  E  ]  +  V  (e 

su  )  /  j  <5  rtuv  \  stuv 

tuv 


(363) 


where  Prs is  the  permutation  operator  that  exchanges  orbital  r  with  orbital  s  .  Although 
the  unitary  generators  (and  hence  density  matrices)  are  not  symmetric,  only  the 
symmetric  portion  of  those  matrices  will  appear  in  the  product  with  the  symmectric 
matrices  hand  g.  Hence, 
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Let  us  define  the  Fock  matrices  [3]  [5] 


f1[z]  sy  h[z]D[z] 

rs  ru  us 

u 

F2[z]  =  y  s[z]d[z] 

rs  o  rtuv14'  stuv 

tuv 

Fiz]  =  Fl-”[z]  ,  fAz] 

rs  rs  1  rs 


(364) 


(365) 


which  allow  the  orbital  gradient  vector  to  be  simplified  as 

f°  =2(f}z]-F}z])  (366) 


Since  an  analytic  derivative  is  never  taken  of  this  vector,  it  is  not  necessary  to  further 
transform  it. 

To  determine  the  elements  of  the  z(R)  vector,  Z  (/?),  we  note  the 
transformation 
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(367) 


Dff(«)  =  [“P(-z(R))DW(R)“P(z(R))]„( 

Using  the  Baker-Campbell-Hausdorff  theorem  again,  we  expand  the  exponentials  and 
find  the  relationship 

<36s> 

where  higher-order  commutators  are  ignored.  We  can  break  this  commutator  into  its 
two  component  terms,  which  we  can  then  express  as  summations: 

D[S(R)  =  ^(RVI.{D[S(R)Z„(R)-Zpr(R)D^\R))  (369) 

Let  us  take  the  gradient  of  this  equation  and  evaluate  it  at  the  reference  geometry: 

d[2  w = <  k)" + z(<]  (*o)(z„,  («„)')-(z„  («„rK'  w)  otoi 

where  we  have  used  the  fact  that  Z(i?0)  =  0 .  To  complete  the  evaluation  of  this 

equation,  consider  that,  within  the  invariant  subspace  of  interest,  the  basis  was  resolved 
by  diagonalizing  the  density  matrix  [3].  Then  the  equation 

DJf[z]  =  DJ'[Z]S  (371) 

pq  pp  pq  V  / 

must  hold.  When  this  equality  is  recognized,  equation  (370)  becomes 

Evaluation  for  /z^gyields 
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(373) 


0 = dW  w +«'  w- Dl‘J(R„))z„(R„y 


z„(K)‘  = 


D^(R0)-D^(R„) 


(For  p  -q,  equation  (372)  is  not  useful  for  identifying  Zpq  (R0)x  )■ 

Diagonalizing  the  density  matrix  is  not  the  only  method  for  resolving  bases.  The 
diagonalization  of  Fock  matrices,  such  as  those  defined  in  equation  (365)  as  well  as  the 
Fock  matrix  defined  by 


eg1  (R)  =  ( R) + I[(2*S,  (*)-  *2,  ( R))  off1  (R) 

rs 


(374) 


may  also  be  used.  The  analogue  of  equation  (367)  is  true  for  these  Fock  matrices  as 
well.  Thus,  depending  on  the  method  of  resolution  used,  we  find  that 


Z„(*o) 


x 


(375) 


and 


Z„(*o) 


x 


gig  w 

e;g(R„)-fiig(R„) 


(376) 


as  well.  Each  of  these  may  contribute  to  the  DCT,  and  so  we  find  that  the  vector 
product  z(Rq)'  has  contributions  equal  to 

(?  w  •  7s  (k))° = ?  w  <  (377) 

(4^)'  -/,2  K)f = 2%'W(*T  K  (378) 

p*q 

and 
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(379) 


(4^)*-7,SW),,=Ee!?W^ 

p*q 

where 


(380) 


is  not  differentiated,  and  can  be  calculated  as  needed.  The  first  factor  under  each 
summation,  which  is  a  derivative,  must  be  further  transformed  into  the  atomic  basis. 
The  gradient  of  each  Fock  matrix  will  involve  gradients  of  one-  and  two-electron  integral 

matrices  (h',gc)  and  gradients  of  density  matrices  (D*,d*)  (see  equations  (374)  and 


(365)),  while  transforming  back  to  the  [5]  basis  will  involve  a  gradient  of  the  K (i?) 

matrix  for  both  Fock  and  density  matrix  gradients.  This  transformation  and 
differentiation  will  yield  three  types  of  terms:  those  involving  gradients  of  integral 

matrices,  involving  gradients  of  K (i?),  and  involving  gradients  of  density  matrices. 

The  terms  involving  gradients  of  the  integral  matrices  will  take  the  form 

Tr  (h^  (R0  )x  Dx  )  +  ± Tr  (g^  (R0 ) '  dx  )  (381) 


where 


K)- a«zjM(r0)) 
DF  -i{D^(7?o);AF} 
d'=ijd|!1(^;A') 
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in  which  the  operation  {;}  is  the  matrix  anticommutator.  The  terms  involving  gradients 
of  the  K(i?)matrix  will  all  taketheform 


k(Roy-fax  (383) 

which  can  be  added  ontothesimilarterm  in  equation  (360).  The  terms  involving 
gradients  of  density  matrices  will  take  the  form 

(384) 

where  p(R)  is  the  vector  of  parameters  for  the  rotation  of  CSF  coefficients  (see 
equation  (167))  and  fcsf  is  analogous  to  fCI ,  but  where  the  commutator  is  taken  with 
the  generators  of  exp(p) .  Now  the  vector  product 

(  fCI  _i_  fCI Q  i  faF  ' 

(*  W'  P(^')')[7c,0  +7c,pJ  '385' 

can  be  cast  as  [3] 

Tr  (h^  (R0 )x  (Da  +  Dk  ))  +  j Tr  (g^  (R0 )*  (dA  +  dk  ))  (386) 

Thus  all  terms  in  equation  (360)  are  written  as  traces  of  gradients  of  integral  matrices 
multiplied  by  various  density  matrices.  Substituting  the  results  of  (381)  and  (386)  into 
equation  (362)yields 

(r!Z|  (Re )  | ff 121  («0 )'  |  rlz|  (R, ))  =  Tr  (h|sl  ( R„ )’  ( D[s|  ( R, )  +  Dq  +  DF  +  DA  +  Dl )) 

+i7/(gl',(«jldM(ffll)1d°  dF+d'+d‘))  (387) 

=  7>(h|sl  (R„y  D[fl  (/?„))  +  4rr(g[s|  («„)'  d£  («„)) 
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These  integral  matrices  have  analytic  forms  in  the  original  atomic  basis  to  which  we 
should  like  to  transform  them.  Recall  that  the  transformation  between  the  [C]  and 
[ S ]  bases  is 

h[s]  ( R )  =  S[C]^  (R) h[c]  (R) S[C]~^  (/?)  (388) 

It  follows  then  that  the  derivative  of  this  equation  evaluated  at  the  reference  geometry 
is 


h[S|  (R)'  =  ff™  (R)*  h|c|  (R)SF™  (R) 


,[c] 


[c]. 


,[c]- 


+  S[c]-/2  (R0)h[ci  (Rj  SlLT/2  (R0)  +  SLL'r/2  (tf0)h[CJ  (R0) s[tr/2  ( R0 ) 


.Mi 


,[c] 


[c]- 


.M  i 


(389) 


Using  equation  (161),  we  find  that 


S[c]"^(/?0)  =  1 

s[c^K)‘  =  -is[c](^)- 


(390) 


and  equation  (389)  is  simplified  as 

1>[S1  (*.)'=-  iS[cl  (R  Y  h[c>  (R )  + 1 hM  (R )'  -  {  h[cl  ( R )  S[cl  ( R )' 

(391 

=  hlcl(R,,)’-i(slcl(R0)',h[cl(R0)} 

Similarly, 

gW  (R)“  =  gM  (R)’  -i{sM  (R)'  .g[<r|  (R)}  (392) 

leading  to  the  transformation  of  equation  (387): 
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(393) 


¥f](R())  H[zl(Roy  ¥f\R0))  =  Tr  h[c\R0)xV[yj(R0) 


rM, 


Jcli 


+  irr 


--Tr 
2 1 ' 


«|C|(*o)'CC|(«.) 


y/[c] 


h[c](/ii),Slcl(«0)’|-D1',(^) 


j[c] 


»lcl, 


-iyy 

4  -1  ' 


( if,, ) ,  S LC|  ( if,, ) ' } .  d£'  ( R„ ) 


:[c] 


i  Mi 


in  which  we  have  used  the  fact  that  the  density  matrices  in  the  two  bases  are  the  same 
at  the  reference  geometry.  Let  us  simplify  the  arguments  of  the  trace.  We  have,  using 
simplified  notation: 

7>({h.S}  D)  =  7r(hSD+ShD)  =  V|,  s  l»  S  I,  l>  (394) 

ijk 


Since  hand  S are  symmetric,  this  is 

I'l  s  n.  Sh  n  =  vs  I,  »  +ssh#D,,  (395) 

ijk  ijk 

The  density  matrix  is  not  symmetric;  however,  since  h  and  S  are,  only  the  symmetric 
portion  of  Dcontributes  to  the  trace.  Thus 

Xs  h  »  ■>  h  ■>  +Ssh„D„  =22>(SF')  (396) 

ijk  ijk 

where  the  Fock  matrices  were  defined  in  equation  (365).  If  we  apply  this  derivation  to 
the  terms  in  equation  (393),  we  find  that 


1 

2 


Tr 


Ac] 


-Tr 


S[C]  W'F, 


1[C] 


(397) 


and,  through  a  similar  derivation, 

Equation  (393)  can  consequently  be  written 


-iy> 

4  ' 


^(RQ),S^(RQy\.^(RQ) 


:(c] 


iMi 


=  -Tr 


(398) 
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(399) 


{vf]  (Ro)|«|z|  (*„)' \vf]  (Re)) =7>[hlcl  («„)*  Dg  (R0) 

+!7v[g[cl  (SoJ'dE1  («0)]-7v[s|c|  (R,y  FLC1  (R„)_ 

The  transfer  to  the  atomic  basis  is  relatively  simple.  Any  operator  in  the  [C]  basis  can 

be  expressed  in  the  [ %\  basis  by  the  transformation 

A[c]  (R)  =  CT  (Rq)  AW  {R)C(R0)  (400) 

(see  equation  (157)).  Since  C(R{))  is  constant  with  respect  to  nuclear  coordinates, 

A[c]  (. R)x  =  CT  (Rq)  AW  (R)x  C (Rq)  (401) 

however,  density  matrices  transform  as  [3] 

D[c]  (R)  =  CT1  (R0)Dm  (R)( C  1  (/?0))r  (402) 

Thus  equation  (399)  can  be  restated  as 

(r!Z)  ( «« )  |  ■ [Zl  («„ )'  |  r!Zl  («0 ))  =  Tr  [cr  ( R0 )  hM  («„ )'  D«  («0 )  (C-1  )Z  (R0 ) 

+  ¥■ Tr  [cr  (K )  SW  ( K  T  (*o )  (C-‘  f  (R0  )]  (403 ) 

-  Tr  [cr  ( R„ ) Sw  ( Ra )'  F«  («„ ) (C"‘  f  (R, ) 

Since 

7r(AB)  =  7r(BA)  (404) 

it  follows  that  the  last  matrix  in  each  trace  argument  can  be  brought  to  the  front,  and  so 
the  transformation  has  no  effect  on  the  form  of  the  trace: 

(r‘z|(^)|»|zl(^)'|r!z|(^))=7'r[(fl0)hW(i;0)'DW(fl0)‘ 

_  r  -  (405) 

+  j7>  gUi(K)‘i2(K)  ~Tr  SW(R0)'F,W(R0) 
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Thus  the  gradient  is  now  written  with  derivatives  only  of  matrices  whose  elements  are 
known  analytically,  i.e.  derivatives  of 

bf(R)  =  \drx'(r,R)  Zj(nR) 

a  I  Ka  r\_ 

(*)=//  drAi  Xi  (r{,R)x*k{r2\R)  r 

Uri  r2\ 

sW  (R)  =  \drx’(r,R)x,(r.R) 

while  the  density  matrices  can  be  calculated  as  needed.  Thus  there  is  no  need  to  use 
finite  differences  to  estimate  the  gradient,  nor  to  create  analytic  functions  of  the  Cl 
coefficients. 

The  Cl  DCT  is  derived  in  the  same  manner,  with  the  exception  that  transition 
density  matrices  are  used,  and  the  entire  quantity  must  afterward  be  divided  by  the 
energy  difference: 

AE K )  fS  (R« y  =  (t^1 ( R, ) I H[zl  (K )“  v?1  ( R„ )) 

=  J>[(^)hW(*o)'D^1(«0)]  (407) 

+ *7>  [gw  (R,  y  («,,)]-  Tr  [sw  (R, )’  F,'«  (R,)_ 


Xj(ri>R)Zi(r: 2,R)  (406) 
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Appendix  H.  The  CSF  DCT 


The  CSF  DCT  is  [5]  [49] 

fT  (*)'  ^  ( R)c‘, '  (R){<Pj  (r :  tf)| (R)}  (408) 

ij 

(see  equation  (244)  in  the  main  text).  Unlike  the  Cl  DCT,  in  this  case  the  gradient  is 
taken  of  the  basis  function.  Since  the  Gel'fand-Tsetlin  basis  functions  can  be  written  as 
linear  combinations  of  Slater  determinants,  let  us  consider  the  gradient  of  a 
determinant: 

it  I  MbM)  =  |  it  it  )  Me  -}  +  |  <t>a  ~k  k  ")  +  |  Mb  A  (  &  )-)  +  - ■  (409  ) 

The  n  !n  terms  in  this  expansion  each  contain  nuclear  derivatives  of  functions  of  all 
electronic  variables;  however,  one  can  also  separate  the  terms  into  a  set  containing 
nuclear  derivatives  of  functions  only  of  electron  1,  a  set  containing  nuclear  derivatives 
of  functions  only  of  electron  2,  etc.  For  example,  consider  the  simple  Slater  determinant 

Its  nuclear  gradient  can  be  written  as 

(411) 

= (1)) A  (2))“l*(^  W)^-  (2))+ k-  (!)*(^  (2)))-|  A  W^r(A  (2))» 

which  can  be  separated  into  two  terms  as  discussed  above.  Thus  we  can  recast  the 
nuclear  derivative  as 
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_8_ 

8R. 


(412) 


i 

where 

*(0|A  =  |  afc&  (?';^))^  (413) 

Since  -J-  is  then  effectively  a  one-electron  operator,  it  can  be  represented  in  the 

molecular  orbital  basis  by  the  product  of  generator  matrices  with  integral  matrices 
much  like  the  one-electron  Hamiltonian  (see  equation  (31)) 

(414> 

r,s 

where 

(415) 

We  can  substitute  this  definition  into  definition  (408): 


^M))-S4(*b(«X«(*)IZ^(*(*#(0K(':«))k(*f 


= X 14  (*)  W*)l  (*)  k,  W)  {*.  Ml  s\i  (k)) 

rs  ij 

=I(^M|4|^M)(a(«)|%M) 

rs 


=  1  o"[z'(«)(#1(«)  *#>(*)' 


(416) 

where  all  pieces  can  be  assumed  at  this  point  to  be  expressed  in  the  [z]  basis.  Thus, 
rather  than  acting  on  the  Gel'fand-Tsetlin  basis  functions,  we  can  assume  the  gradient 
acts  upon  the  molecular  orbitals  directly.  We  can  now  transform  this  form  of  the  CSF 
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DCT  back  into  the  [C]  basis  just  as  we  have  done  with  the  Cl  DCT.  The  derivative  of  a 
MO  evaluated  at  the  reference  geometry  can  be  expressed  as  the  transformation 


dR, 


d R, 


(417) 


#(*»))  =  £  exp(K(«))exp(Z(«))|#l(«)) 

'(l  +  K(*)+...)(l  +  Z(*)+...)|*!I1(tf))’ 

jl  +  K(fl)+Z(tf)  +  ..,)|<if  !(*))' 

where  all  other  terms  contain  at  least  one  factor  of  K  (/?n)orZ(/?0)  and  thus  evaluate 
to  zero.  Transforming  from  the  [S]  basis  back  to  the  [C]  basis  is  effected  by 


dR, 


(*0  ))  =  ST  (s[Cl_K  («o  ))  |  i 1  (*»))  +  SW"X  («„  )  A I  ff '  («o  )) 


(K  ( K  ))S[r^  (R„)  |  # 1  ( )}  +  £  ( Z  (*,,))  S[r|^  (R0 )  |  («0 )) 


8R, 


(418) 


(see  equation  (163)),  which,  using  equation  (390),  simplifies  to 


dR, 


«f 1 K ))  =  -  i  w,  (s[cl  K ))  I  rfc|  (K))+ w,  I  «f'  (*o )) 

+7lr(K(«.))|#1(«»))+J:(z(^))|#l(^)) 


(419) 


Since 


X], 


(«»)  ={#'(«»)  =(#'(fio)  =wqw 


pi, 


,[C], 


(420) 


it  follows  that  integrating  equation  (419)  against  (^ZJ(/?0)|  yields 

+(#'(«.)|i(K(R.))|#1K))+(«iJc|(«»)|i(z('*o))|#1(«.)) 


(421) 


Elements  of  ^-S^(/?0)are 
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( A' S|cl  («0  ))„  =  {&  (f 1  (Ro)  I  «f 1  («0 )}  ■ + (tf1  K )  I  i  «f 1  ( K )} 


(422) 


and  so  the  first  two  terms  on  the  right-hand  side  of  equation  (421)  can  be  combined: 


#’(«.)  a:<f1K)H(WclK)  *  4C1  («o ))  -  (Jr  4C1  ( )  <f]K) 


+ 


#'(«.)  +(K(«„))  «f'M+((f'(*o)  *(z(fl.))  #'(«») 


(423) 


If  we  use  the  parameter-generator  form  of  K  and  Z  (see  equation  (355)), 

KK)=XW*o)(E,v-£,v) 

r'<s' 

ZW=Ev,'ft)(J,v-«,v) 

r'cs' 


(424) 


then  the  final  terms  in  equation  (423)  become 


am, 


am, 


#]K)  *(ZK))  #](^»))  =  Z^(^vK))(#]K)  (Er.,-E,r)  #](*0) 


aM, 


aM, 


(425) 


as  only  the  parameters  (£ryand  zrV)  have  nuclear  dependence,  while  the  generators  ( 


ErW  and  EsV )  participate  in  the  integral.  Because  of  the  annihilation/creation  form  of 
the  generators  (see  equation  (30)),  those  integrals  become  delta  functions: 


#]K)  ^(k(*0))  rfc]K))  =  S4r(^vK))(^A,-^Av) 


am, 


4C]  (Ro)  MZ(Ro))  #]  K)}  =  ZMzr,  (R0  ))(Sr.rS,'-8r.'8,r) 


am. 


(426) 


Because  of  the  restricted  sum,  only  one  of  those  delta  products  will  ever  be  non -zero; 
thus  these  equations  reduce  to 


(f'K)  J7(kK))#]K) 


,[c] , 


#]  W  MZM)  #]K) 


w, 


dRx  krs  (  R()  ) 


(427) 
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The  various  forms  of  ^§- zrs(RQ )  were  discussed  in  equations  (373)-(376).  From 
equation  (423)  let  us  define 

(which  is  antisymmetric)  so  that  the  entire  CSF  DCT  becomes 


(428) 


X  c*  (Ro )  cjj  (Ro )  \<Pt  ( r;  rM  v,  (r;  Ro ))  = 

ij 

ZD"[C|  K)(/„*T,[cl K)' +wKAK)+ («„)) 

r<s 

=  Tr ( D"[c>”  ( R„ )  fZm  (R„ ) ' )  +  Z  DT'  (*. )  £ K  (R« )  (429) 

r^s 

+Zi:B"w(«»)Ar'JD+Zi^w(^)A?'"F 


+Zie:f  KM 


CSF  IJ  Q 


where  we  have  defined 


aCSFIJX 

pq 


dT(r o) 


(430) 


(see  equation  (380))  and  the  antisymmetric  part  of  the  one-electron  density  matrix  is 
sufficient  for  the  trace  with  the  CSF  orbital  gradient  vector.  The  last  three  terms  are 
treated  analogously  to  the  terms  in  equations  (377)-(379),  with  the  appropriate 

substitution  of  the  A^f  //xterm  The  second  term  joins  them  in  a  manner  similar  to 


(385),  with  the  resulting  form 


jyjhM  (/;„)' (» 

I7v(gl<:|(^)'(d 


+  2 


CSFIJ F  pCSf  IJQ  Y)CSFU  A  _|_  jjCSFWJi. 

CSFIJ F  _|_(jCSF/./Q  +  ^CSFIJA  +  ^CSFIJX 


(431) 
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where  DCSF/yxand  dc?/  //x  are  defined  analogously  to  those  found  in  equation  (387). 
These  traces  can  now  be  transformed  into  the  atomic  basis  similar  to  equation  (403): 

Z  c/'  ( R<> )  cJj  (Ro){^{r>Ro)\w;iPj{r>Ro))  = 

ij 

rr(D'yW"(^)/“FU(^)') 

+7>(hw  («„)’  (D™  +Dcs'"q  +Dof“a  +  Dcsf"1)) 
+I7'r(gw (R„f  (d“F"F  +dof"5  +dCSF"A  +d°F"1)) 

=  rr(D"w"(fiD)/“FU(^)‘) 

+  rr(hW(ff„y  D“™)  +  iTr(gW  (R„)‘  dFf") 

At  this  point,  the  CSF  DCT  has  been  defined  in  terms  of  density  matrices  and 
gradients  of  integral  matrices  at  the  reference  geometry.  Combining  this  with  the  Cl 
DCT  term  gives  the  whole  DCT  in  terms  of  functions  which  can  be  differentiated 
analytically. 
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Appendix  I.  Angular  Momentum 


Definition 

Classically,  we  consider  angular  momentum  as  the  cross  product  of  an  object's 
linear  momentum  with  its  position  with  respect  to  some  fixed  point  [25]: 

J  =  rxp  (433) 

Since  it  is  the  cross  product  of  two  polar  vectors,  angular  momentum  is  a  pseudovector; 
that  is,  it  rotates  like  a  polar  vector,  but  it  remains  invariant  under  inversion,  unlike  a 
polar  vector.  However,  in  this  paper  we  will  be  dealing  with  special  unitary  groups 
which  do  not  allow  inversion,  so  angular  momentum  can  be  treated  summarily  as  if  it 
were  a  polar  vector. 

As  with  the  Hamiltonian,  we  will  derive  the  quantum  mechanical  equivalent  of 
equation  (433)  by  turning  the  classical  quantities  into  operators  [61]: 

J  =  r  x  p  (434) 

Recalling  the  quantum  mechanical  definition  of  the  linear  momentum  operator,  (284), 
we  see  that  the  angular  momentum  operator  is,  component -wise  [46], 

J  (435) 

The  formulation  in  equation  (434)  was  used  in  Appendix  B  when  exploring  the  spin-orbit 
Hamiltonian,  where  L  was  the  orbital  angular  momentum  of  the  electron.  However, 
there  we  also  found  the  set  of  spin  operators  &i  of  an  electron  was  considered  a  set  of 
angular  momentum  operators.  An  electron,  though  massive,  is  considered  a  point 
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particle  in  quantum  chemistry,  so  it  is  a  fallacy  to  attempt  to  connect  this  spin  angular 


momentum  with  that  of  a  spinning  orb,  and  equation  (434)  does  not  apply.  What  we 
need,  then,  is  a  more  general  definition  of  angular  momentum. 

In  an  intuitive  sense,  angular  momentum  generates  rotation.  In  a  mathematical 

sense,  angular  momentum  operators  are  the  generators  of  rotation.  A  rotation  U  in 
three-space  can  be  expressed  as  the  exponentiation  of  it  rotational  parameters 
multiplied  by  the  generators  [34]: 


U  =  exp(^r/v  +  0yJ  y  +  6ZJZ )  (436) 

This  is  certainly  not  the  only  parameterization  of  a  three-dimensional  rotation, 
but  it  is  as  valid  as  any.  The  set  of  rotations  is  a  Lie  group,  and  so  the  set  of  angular 
momentum  operators  forms  its  underlying  Lie  algebra.  Lie  algebras  are  defined  by  their 
commutation  relations  [34],  In  this  case,  the  operators  obey  the  relations 


J  J 

x  ’  y 


Jy’Jz 


J  J 

z  ’  x 


=  iJz 

=iK 

=  i  J , 


(437) 


which  is  a  specific  case  of  the  general  Lie  algebra  commutation  relation 

A.OJ]  =  2>wd,  (438) 

i,j,k 

It  is  easy  to  see  that  the  angular  momentum  as  formulated  in  equation  (434)  obeys  the 
relations  in  equation  (437)  given  the  fundamental  commutation  relation  between 
position  and  linear  momentum, 
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[x,  px  ]  =  /,  etc. 


(439) 


however,  we  say  that  any  set  of  operators  that  obeys  the  commutation  relations 
presented  in  equation  (437)  are  angular  momentum  operators  [46],  Thus,  even  though 
the  spin  operators  do  not  conform  to  the  formulation  in  equation  (434),  they  do  obey 
the  proper  commutation  relations  and  thus  can  be  considered  angular  momentum 
operators. 

Eigenvalue  equations 

As  angular  momentum  is  an  observable  quantity,  it  obeys  an  eigenvalue 
equation.  In  this  section  we  will  discover  the  possible  eigenvalues  and  eigenfunctions 
associated  with  the  angular  momentum  operators. 

Commuting  Operators 

It  is  clearfrom  equation  (437)  that  none  of  the  angular  momentum  operators 
commute;  this  implies  that  only  one  of  the  quantities,  Jx,  Jy ,  or  J . ,  can  be  observed 

precisely  at  any  given  moment  of  measurement  [26],  This  also  implies  that  they  form  a 
rank  1  algebra  (i.e.,  the  largest  abelian  subalgebra  is  one-dimensional)  and  thus  has  one 
Casimir  invariant  operator  [34]: 

J]  +  J]  +  J)  =  J2  (440) 

(This  operator  is  discovered  by  solving  the  secular  equation  for  an  arbitrary  element  in 
the  algebra;  this  process  is  discussed  in  Gilmore  [34]  but  I  will  not  elaborate  here.)  The 
Casimir  invariant  is  a  multiple  of  the  identity  matrix,  and  thus  commutes  with  all 
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elements  of  the  algebra.  This  further  implies  that  any  one  of  the  aforementioned 
quantities  can  be  observed  concurrently  with  the  Casimir  invariant;  typically  we  choose 
J  _ .  Thus  we  have  a  maximal  set  of  two  commuting  observables  in  the  set  of  angular 

momenta,  J_  and  J2,  which  share  a  set  of  eigenvectors.  Geometrically,  this  implies 
that,  while  we  cannot  know  the  precise  direction  of  the  angular  momentum  vector,  we 
can  know  its  length  and  its  projection  onto  the  z-axis. 


Eigenvalues 

Thus  we  have  the  following  eigenvalue  equations: 

J 2 1  jm)  =  Aj  |  jm) 

Jz  |  jm)  =  4  |  jm) 


(441) 


While  the  operations  7r  |  jm)  and  jy\jm)  are  of  lesser  interest,  the  operations  j+\jm) 


and  J  |  jm) ,  where  we  define 
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(442) 


K  =  K + iK 

J  =l~iJy 

which  are  adjoints  of  each  other,  are  of  particular  use.  Using  these  in  place  of  the  Jx 
and  J  operators,  we  have  a  new  Lie  algebra  with  commutation  relations  [46] 

JZ,J±  =±J± 

~J+J_]  =  2Jz  (443) 

”  A  ~  /V  “ 

J~,J±  =0 

Since /"and  J ±  commute,  the  eigenvalue  equation 

/ 2/±  |  jm)  =  AjJ±  |  jm)  (444) 

implies  that  J±\jmj  is  also  an  eigenvector  of  J2  with  eigenvalue  /L ,  while  the  equation 
JJ±  |  jm)  =  (+/±  +  JJz ) |  jm)  =  (4  +  l)/±|  jm)  (445) 

implies  that  it  is  also  an  eigenvector  of  J _  with  eigenvalue  ( Xm  ±l) .  Furthermore, 
equation  (445)  implies  that 

J ±  |  jm)  =  C±\j m ± l)  (446) 

which  shows  that  the  eigenvalues  of  J z  differ  by  integer  steps.  Combining  definition 
(440)  with  definition  (442),  we  find  a  new  form  of  the  Casimir  invariant, 

J2=)jJ  +{J  J++j2z  (447) 
We  can  takethe  expectation  value  of  this  equation, 
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(448) 


(jm  I  J2 1  jmj  =  \\jm  I  jJ  I  jmj  +  \(jm  I  J  J+  I  jrnj  +  (jm  I  /:  I  jtnj 

Aj=?(\C-  |2+|C+|2)  +  ^ 

which  implies  the  quantity  A} must  be  positive.  Since  both  A.  and  Am  are  real, this 
further  implies  that  A.  is  non-negative  and  that  Am  is  bound  above  and  below  by  ±-/l~ 

J  m  y  J 

[69],  Let  jmnax )  be  the  vector  whose  eigenvalue  of  is  maximal.  Then  it  is  true 
that 

K  |Mnax)  =  0  (449) 

and  hence  that 


TT  |  ) = (i2  -  T  -  7  )|  Mj = 0 

Expanding  this  last  equation,  we  find 

+1) 

Thus  we  conclude  that  for  maximal  Am, 


By  analogous  argument,  for  the  minimal  value  7L 

mmin 

A.=A  (A  -l) 

J  '"min  V  '"min  / 

Comparing  equations  (452)  and  (453)  leads  to  the  conclusion  that 

mmdx  mm\n 


(450) 


(451) 


(452) 


(453) 


(454) 
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Since  the  values  of  Am  differ  by  integer  steps  (see  equation  (446)),  this  conclusion 
implies  that  Am  must  take  on  integer  or  half-integer  values.  Ifwelet  Am  be  assigned 

the  symbol  j ,  then  the  eigenvalue  of  J2  is  j(j+ 1)  and  if  we  let  Am  be  assigned  the 
symbol  m ,  m  ranges  from  -  j  to  j  .  The  equations  are 


J2\jm)=  j(j  +  l)\jm) 
Jz\jm)  =  m\jm) 


(455) 


Thus  not  only  is  the  z-projection  quantized,  it  is  limited  by  the  total  angular  momentum 
(see  Figure  44).  Table  4  shows  the  first  several  values  of  j  and  the  possible  values  of  m. 

Apart  from  the  eigenvalue  equations,  we  also  have  the  equations  involving  J+  and  J  . 


Usingthe  relationship 


j2  =  j±L+  +  j:  +  jz 


(456) 


Figure  44.  Angular  momentum  allowed  observables 
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Table  4.  Angular  momentum  values 


j 

m 

0 

0 

1/2 

-1/2 

1/2 

1 

-1 

0 

1 

3/2 

-3/2 

-1/2 

1/2 

3/2 

2 

-2 

-1 

0 

1 

2 

We  find  through  integration  that 


(jm |  J 2 1  J/n)  =  (.//« |  J±iT  |  jm)  +  (jm  \  ]\  \  jm)  +  ( jm  \  J.  \  jm ) 
j(j  +  l)  =  \CT[  +  m(m+l) 


(457) 


which  leads  to  the  equations 

J±  |  j  m)  =  yjj(j  +  l)  —  m  (m±\)\j  m  ±  l)  (458) 

The  eigenvalue  equations  (455)  paired  with  the  equations  immediately  above  will  be 
most  useful  in  forming  the  matrix  representation  of  the  angular  momentum  operators 
in  the  succeeding  section. 

It  is  an  interesting  aside  to  note  that,  although  angular  momentum  may  have 
integer  or  half-integer  values,  the  two  sets  never  coincide.  That  is  to  say,  a  system  with 
integer  spin  can  never  bemadeto  have  half  integer  spin,  and  vice  versa.  Asystem  with 
an  even  number  of  electrons  can  be  a  singlet,  triplet,  quintuplet,  etc.,  with  the  available 
spin  states  only  being  limited  by  the  number  of  electrons;  however,  it  will  never  be  a 
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doublet,  quartet,  sextet,  etc.  The  opposite  is  true  for  odd-electron  systems.  Fermions 


may  be  turned  into  one  another,  but  can  never  be  transformed  into  bosons.  Orbital 
angular  momentum  only  manifests  itself  in  integer  steps  (see  the  spherical  harmonics, 
below). 

Physical  Eigenfunctions 

The  eigenvalue  equations  (441)  form  solvable  differential  equations  when  the 
operators  J2  and  J _  are  derived  from  equation  (435)  and  the  eigenvalues  are  replaced 


with  j(j  + 1)  and  m.  Rather  than  Cartesian  coordinates,  it  is  easiest  to  solve  these 
equations  in  spherical  coordinates.  By  representing  these  equations  in  physical  space, 
we  are  limited  to  representing  angular  momenta  which  reside  in  physical  space  (e.g., 
orbital  angular  momentum).  This  representation  is  limited  to  integer-value  angular 
momenta,  and  so  will  not  include  spin  (which  does  not  reside  in  physical  space).  I  will 
not  belabor  the  transformation  nor  the  solution  at  this  point;  suffice  it  to  say  that  using 
the  relationships 


x=  r  sin  6  cosy? 

y  =  r  sin  6  cosy?  (459) 

z  =  r  cos# 

the  angular  momentum  operators  become  in  spherical  coordinates  [46] 


i  ■  & 
J  — x  — i  — 

d(p 


1 


a2  i 

-+- 


sin ~  6  d(p~  sin  6  86 


e 

do) 


(460) 
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When  these  are  substituted  into  their  respective  eigenvalue  equations,  they  form  the 


familiar  differential  equations  whose  solutions  are  the  spherical  harmonics,  Ylm  (&,<p),  a 
complete  set  of  orthonormal  functions,  whose  formula  I  include  here  for  completeness: 

,9)  2'n 

Because  the  spherical  harmonics  are  complete  and  orthonormal,  they  provide  a  very 
useful  basis  from  which  to  build  wavefunctions.  When  multiplied  by  a  radial  function, 
these  functions  are  able  to  span  thethree-dimensional  wave  functions  of  single 
electrons;  the  direct  product  of  N  sets  of  these  functions  can  form  a  physical  basis  for  a 
multi-electron  system. 

Rotation  of  Eigenfunctions 

When  a  rotation  is  applied  to  a  spin  eigenfunction,  it  becomes  a  linear 
combination  of  other  spin  eigenfunctions.  However,  these  rotations  do  not  allow  mixing 
outside  the  irrep  to  which  the  eigenfunction  belongs;  that  is,  the  total  spin  does  not  mix 
[46]: 

R\jm)  =  YjKm  (^)|  Jm')  (462) 

m 

where  the  D]m,m  [r^  are  known  as  the  Wigner  D  matrices  [37], 


2/  +  1  ( l-m)\ 
2  (l  +  m)\ 


sin'”  6 


d 


d  (cos  6) 


(cos9  6>-l)  exp  ( iincp^ 


(461) 
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Matrix  Representation 


The  commutation  relations  in  equation  (437)  define  (to  within  a  constant)  the  Lie 
algebra  su(  2)  [34];  thus  the  physics  of  angular  momentum  is  inexorably  connected  to 

the  mathematics  of  the  su( 2)  algebra.  We  gain  great  satisfaction  from  this  association 
as  the  physics  of  rotation,  generated  by  angular  momentum,  and  the  group  of  rotations, 
S6>(3),  (one  of  thegroups)  generated  by  su( 2),  should  be  connected,  as  we  intuitively 
expect. 


There  are  an  infinite  number  of  trios  of  matrices  that  obey  the  su(2) 
commutation  rules.  The  defining  representation  is  the  set  of  2x2  matrices 


"0  0 

o 

i 

f-i  0" 

J  °y 

(463) 

which  the  reader  will  find  obey  the  relations  in  equation  (437)  where  i  is  replaced  with 
-2.  An  adjustment  of  phase  leads  to  the  Pauli  spin  matrices  [46]: 
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which  form  another  equally  valid  defining  representation  of  su( 2)  generators. 

[Note  that  the  matrices  in  the  collection  (463)  are  antihermitian  (symmetric  imaginary 
components  and  antisymmetric  real  components)  while  those  in  (464)  are  hermitian.  In 
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general,  the  physicist  prefers  hermitian  operators,  as,  by  postulate,  they  represent 
observable  quantities  and  obey  real  eigenvalue  equations;  however,  exponentiated 
antihermitian  matrices  are  unitary,  which  are  more  appropriate  for  transformations. 
Thus  when  one  wishes  to  study  angular  momentum  per  se,  an  hermitian  representation 
is  fitting;  however,  when  one  wishes  to  use  angular  momentum  as  a  means  to  the  end 
of  rotations,  an  antihermitian  representation  may  be  more  applicable.] 

The  regular  representation  is  the  set  of  3x3  matrices  derived  directly  from  the 
commutation  rules  [34]: 

"0  0  0  VO  0  lVO  -1  (T 

00  -1,0  0  0,1  0  0  (465) 

v0  1  0  J  [-1  0  oj  [o  0  0, 

In  the  previous  section,  we  found  a  good  basis  to  work  in,  the  |  jm )  basis.  Using 

equations  (455)  and  (458),  we  can  form  the  matrix  representations  of  72,/T,/+,  and  J 
in  this  infinite  but  discrete  basis.  Figure  45  Shows  the  first  10x10  block  of  each 
operator,  where  zeros  have  been  excluded,  although  the  block-diagonal  form  has  been 
emphasized.  Jx  and  7v  can  be  determined  using  equation  (442).  Each  block  or  irrep 

labels  a  different  value  of  j  .  For  this  reason,  each  block  can  be  used  separately  as  the 
case  may  require.  For  instance,  the  2x2  blocks  (which  are  the  defining  representation) 
can  be  used  to  generate  rotations  on  two-component  spinors  such  as  spin,  while  the 
3x3  blocks  (which  are  a  transformation  of  the  regular  representation)  generate 
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Figure  45.  Irreps  of  angular  momentum  operators 


Figure  46.  Addition  of  spin 
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rotations  of  rank-1  tensors,  such  as  p-states  and  vectors  in  classical  three-dimensional 


space.  In  this  basis,  both  /2and  J_  are  diagonal,  as  expected.  Notice  that  each  block  of 
J2  is  an  identity  matrix  scaled  by  its  eigenvalue  j{j  + 1).  The  elements  of  J_  are  the 
eigenvalues  m . 


Coupling 

Addition 

In  the  previous  section,  we  identified  the  set  of  simultaneous  eigenfunctions  of 
J2  and  J  _ ,  |  jm ) .  If  a  system  has  two  angular  momenta  associated  with  it,  it  is  possible 

to  form  an  eigenfunction  of  both  sets  of  operators  by  taking  the  direct  product  of 
eigenfunctions  from  either  set: 


j{m\)®\j2m2) 


h  h 

m,  ni-, 


(466) 


This  is  an  eigenfunction  for  J2 ,  J2,  /lT,and  J2z  .  Alternatively,  we  can  add  the  two 
angular  momenta  together  and  form  an  eigenfunction  of  ( J{  +  /2j  ,/,2 ,  and 

J2  [46],  While  these  two  formulations  are  equally  valid  and  differ  only  by  a  unitary 
transformation,  the  latter  has  an  advantage  for  systems  in  which  the  angular  momenta 
are  allowed  to  couple.  Forsuch  systems,  theset  of  quantum  numbers  jl,j2,ml,  and  m2 
are  not  good  quantum  numbers  since  they  are  allowed  to  change.  The  good  quantum 
numbers,  i.e.,  the  eigenvalues  that  will  not  change,  will  be  the  total  angular  momentum 
and  its  z-projection.  Even  in  the  non-coupling  case,  this  latter  formulation  forms  a  more 
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compact  notation,  with  only  two  eigenvalues  per  function  to  consider,  rather  than  2N  . 
What  more,  this  latter  formulation  allows  the  system  to  be  easily  expressed  in  the  irrep 
basis  discussed  in  the  previous  section. 

Example:  adding  electron  spin 

Building  a  basis  for  the  GUGA  depends  heavily  upon  spin  adaptation,  which  is  the 
process  of  transforming  eigenfunctions  into  pure  states  of  total  j  [20],  The  following 
example  of  spin  coupling  illustrates  that  transformation. 

Consider  a  system  of  two  spin-1/2  particles  such  as  electrons.  Considering  only 
the  spin  portion,  each  individual  electron  wave  function  can  have  either  of  the  forms 


\s  =  y2,ms  =  y2)  =  \i/ 

=H«) 

\s=y,ms  =-y)=  - 

We  indicatethe  direct  product  merely  by  juxtaposition,  so  that  the  uncoupled  functions 
are 


|  a  a) 
|  afi) 
|  Pa) 

m 

where 


(468) 


\rr')=\r)®\f)  (469) 

Now  consider  transforming  to  spin-coupled  eigenfunctions.  Since  each  electron  has 
spin-1/2,  the  total  spin  can  either  be  l/2+l/2=l  or  l/2-l/2=0  (See  Figure  46).  The  spin-1 
eigenfunction  has  three  possible  z -project ions,  -1,  0,  or  1,  and  thus  forms  a  triplet.  The 
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spin-0  eigenfunction  has  no  z-projection,  and  is  a  singlet.  We  can  indicate  these 
eigenfunctions  as 


11) 

10) 

1-1) 

00) 


(470) 


We  now  need  a  way  to  define  these  spin-coupled  functions  in  terms  of  the  original, 
uncoupled  functions  which  are  more  intuitive.  The  transformation 

I  7m)  =  ^  |  mim2  )  (mim2  |  jm)  (471) 

ml  ,iti2 


has  coefficients  we  define  as 


(n\m2\jm)  =  cjmnhnh  (472) 

known  as  Clebsch-Gordan  coefficients  [46],  which  are  tabulated  in  various  locations. 
Figure  47  shows  the  Clebsch-Gordan  coefficients  for  small  angular  momenta  [70],  Using 
these  coefficients,  we  find  the  following  relationships: 

"/  =  rm)  l10)  =  i(l«/8>  +  l^«)) 

ii-i)=i m  i°o>=i(i«/?>-i/?“» 

This  result  was  gleaned  from  the  l/2xl/2  block  of  Figure  47.  Now  consider  adding  a 
third  electron.  Since  the  first  two  are  already  coupled,  we  are  coupling  an  angular 
momentum  of  1  or  0  to  an  angular  momentum  of  1/2.  Thus  we  use  the  lxl/2  block  to 
find  the  first  set  of  functions,  and  the  second  set  (0x1/2)  is  trivial.  We  can  now  have 
total  spins  of  3/2  or  1/2  (see  Figure  48).  Note  that  there  are  two  ways  to  have  a  spin- 
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1/2  system;  either  coupling  1  to  -1/2  or  coupling  0  to  1/2.  Using  Figure  47  again,  we  find 


the  following  eigenfunctions: 

||f  )  =  \aaa) 

Iff)  =  |j(|aa/?)  +  |a/?a}  +  |/?aa)) 

|ff)  =  i(|  <z/)fi)  +  \/)afi)  +  \/)fia)) 

|{f)=|  ppp) 

r  (47^ 

,ii,_  -%(2\aaP)-(\aPa)+\Paa))) 

I  ^(\aPa)-\Paa)) 

22  [ 

Now  we  have  a  quartet  and  two  degenerate  doublets  which  give  us  the  23  =  8 
eigenfunctions.  This  process  can  continue  indefinitely,  coupling  the  spins  of  as  many 
electrons  as  are  in  the  system  in  question.  As  the  number  of  electrons  increases,  so 
does  the  degeneracy  of  the  multiplets.  Figure  49  shows  these  degeneracies  as  the  spin 
functions  are  constructed  genealogically  [39]  [71].  Each  path  through  the  graph  from 
left  to  right  represents  a  different  spin  function. 

Coupled  rotated  spin  eigenfunctions 

When  we  take  the  direct  product  of  rotated  eigenfunctions  we  can  apply 
equation  (466)  to  equation  (462): 
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Note:  A  square-root  sign  is  to  be  unrlersli 

1+1/2 +1/2  1  0  0 
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Figure  48.  Addition  of  spin  for  three  electrons 


Figure  49.  Genealogical  construction  of  spin  functions 
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(475) 


where  7  ranges  from  k  +  j  to  |fc  -  j'|  and  M  ranges  appropriately  for  each  J  ;  the  final 


term  in  the  last  equation  is  the  Clebsch-Gordan  coefficient. 


Wigner-Eckart  Theorem 

*  J  /V 

For  a  given  value  of  j,  the  set  of  simultaneous  eigenvectors  of  J  and  J_, 
\j~j)’\j~(j -1)),—  \j  j  -l),  \j  j),  form  a  tensor  of  rank  j.  For  instance,  for  y=0,the 
only  eigenvector  is  |00)  and  forms  a  rank-0  tensor  by  itself.  For  7  =  1,  the  three 
eigenvectors  1 1  —  l) ,  |10),  and  |ll)  form  a  vector,  or  rank-1  tensor.  Any  tensor  whose 
elements  T(k,q),  when  rotated,  mix  only  among  themselves  is  known  as  a  spherical 
tensor  operator  [46],  Clearly  this  is  the  case  for  the  tensors  of  spin  eigenfunctions  (see 
equation  (462)).  Each  element  of  the  spherical  tensor  operator,  when  applied  to  a  spin 
eigenfunction,  will  thus  act  like  the  direct  product  of  two  angular  momentum 
eigenfunctions: 


T(k,q')\dJM 


k 

q' 


(476) 
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where  d  may  be  a  collection  of  non-angular  momentum  eigenvalues  which  further 


identify  the  function.  By  completeness  we  can  write  an  arbitrary  eigenfunction 

\d'J'M")as 


k 

J\ 

Jk 

J 

q' 

M'l 

V 

M ' 

k 

J\ 

Jk 

J 

q' 

M'/ 

V 

M ' 

,  d'J'M"j 
t  J'M”)®\d' 


(477) 


Let  us  now  integrate  this  equation  against  the  eigenfunction  |  dKQ) 


{dKQ\d'J'M”)  =  {d\®'£(KQ*i  ^ ^  ^J'M'^j<S>\d')  (478) 


multiply  both  sides  by 
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and  sum  over  all  J  'and  M  " : 
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k  J  \  k  J 

J'M ")(  J'M " 

q'  M  '  \q  M 


(480) 


The  last  integral,  a  Clebsch-Gordan  coefficient,  is  real  and  is  its  own  conjugate: 


UdKQ\dTM^q  JM\rM ")  = 
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Through  completeness,  this  forms  a  delta  function  on  the  right  side: 


£  ( dKQ\d'J'M " 
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while  on  the  left  side,  the  integral  over  the  spin  functions  also  forms  a  delta  function: 


(dKQ\d'KQ) 


Ik  J 

\q  m 


KQ)  =  (d\®(KQ 
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Let  us  substitute  equation  (476)  into  equation(483): 


(dKQ\d'KQ) 


Ik  J 
\q  M 


KQ)  =  {dKQ\T(k,q)\d'JM 


(484) 


and  define 


(dKQ\d'KQ)  =  (dKWTk  II d'K) 


(485) 


to  be  the  reduced  matrix  element  of  the  spherical  tensor  operator  Tk .  Thus  we  have 
that  the  matrix  element  of  a  spherical  tensor. 


(dKQ\T(k,q)\d'JM) 
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(486) 


has  been  separated  into  geometrical  (the  Clebsch-Gordan  coefficient)  and  dynamical 
(the  reduced  matrix  element)  components  [46],  This  separation  allowed  the  reduced 
matrix  element  to  be  treated  separately,  which  led  to  the  definition  of  spin -orbit 
generators  in  terms  of  spin-free  unitary  generators  [44], 
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Appendix  J.  The  Symmetric  Group 


The  symmetric  group  is  the  group  of  all  (finite)  permutations.  It  contains,  as 
subgroups,  the  permutations  of  N  objects,  SN,  each  with  N\  elements.  It  bears 
importance  not  only  because  of  the  permutation  of  electrons,  but  also  beca  use  of 
Cayley's  theorem,  which  states  that  any  finite  group  is  isomorphic  to  a  symmetric  group 
[72], 

Cyclic  notation 

Elements  of  this  group  are  best  expressed  in  cyclic  notation  [72];  for  example, 
the  symbol(l  3  2)  indicates  that  the  first  object  moves  to  where  the  second  object  was, 

the  third  object  moves  to  where  the  first  object  was,  and  the  second  object  moves  to 
where  the  third  object  was  (this  is  the  right-to-left  convention;  a  left-to-right  convention 
is  equally  valid).  A  cycle  containing  only  two  numbers  is  a  transposition ;  a  transposition 
of  the  form  (/  z  +  l)  is  called  an  elementary  transposition.  A  product  of  cycles,  e.g. 

(l  2)(l  3),  results  in  the  permutation  that  takes  1  to  3,  3  (to  1)  to  2,  and  2  to  1.  A  few 

important  theorems  about  the  symmetric  group  in  cyclic  notation  are: 

1.  Any  permutation  can  be  expressed  as  a  product  of  transpositions  or  elementary 
transpositions;  if  the  number  of  transpositions  is  even,  the  parity  of  the 
permutation  is  even  (likewise  for  odd) 

2.  The  set  of  even  permutations  in  a  group  forms  a  subgroup  called  the  alternating 
group 
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3.  Cycles  do  not  generally  commute  unless  they  are  disjoint;  i.e.,  they  have  no 
numbers  in  common 

4.  Any  permutation  can  be  expressed  as  a  non-unique  product  of  disjoint  cycles 
A  number  which  is  not  permuted  may  be  noted  as  a  cycle  by  itself,  or  simply  left  off 
entirely.  For  example,  thetranspositions  (l  2),  (l  2)(3),  and  (l  2)(3)(4)  are  all 

equivalent.  In  this  way,  it  is  clear  that  the  transposition  (l  2)  is  not  only  an  element  of 
S2,  but  also  of  S3,  S4,  ad  infinitum.  In  fact,  every  group  SN  is  a  subgroup  of  SN+l,  per 
the  subgroup  chain 

51<=52cz53cz...cz5jv_1c^  (487) 


Young  tableaux 

The  regular  representation  of  a  permutation  pk  in  the  group  SN  is  an  n\xn\ 
matrix.  The  basis  vectors  for  this  representation  are  constructed  as 


M  =  Pk\l) 

{P*\  =  (1\P? 


(488) 


(the  latter  property  being  true  because  symmetry  operators  are  unitary  [37]).  We  then 
construct  the  elements  as 


(R(Pk)\=(Pi'P^Pj) 


1  if  P?PkPj=l 
0  otherwise 


(489) 


In  this  form,  the  identity  element  will  be  diagonal,  but  in  general  the  other  permutations 
will  have  no  structure.  The  irreducible  representation  basis  vectors  are  specific  linear 
combinations  of  the  vectors  in  definition  (488)  which  can  be  labeled  by  Young  tableaux 
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[39],  Each  Young  frame  of  N  boxes  labels  an  irrep  of  SN .  For  the  group  SN,  we  take 


the  N  objects  to  be  permuted  (usually  the  numbers  1...N)  and  disperse  them  into  each 
of  the  Young  frames  of  n  boxes  using  the  following  rules  [39]: 

1.  No  number  is  to  be  repeated 

2.  The  number  must  be  largerthan  the  number  to  the  left  or  above  it 

The  frames  filled  with  numbers  by  these  rules  are  standard  Young  tableaux,  and  each 
one  is  a  basis  vector  for  the  irreducible  representation.  Figure  50  explicitly  shows  all  the 
Young  tableaux  that  label  S5  irrep  basis  vectors.  Each  basis  vector  is  labeled  by  a  pair  of 
those  tableaux,  but  the  pair  must  havethe  same  frame.  Thus,  if  k  tableaux  have  the 
same  frame  in  a  group,  they  will  be  used  to  label  k 2  basis  vectors.  Figure  51  shows  an 
example  basis  vector  label  using  two  tableaux  of  the  same  shape.  Using  Figure  50,  we 
see  then  thatthe  S5  irrep  will  havel2  +  42  +52  +62  +52  +42  +12  =120  =  5!  basis  vectors, 
the  same  as  the  regular  representation.  In  the  SGA,  each  box  represents  a  spin-orbital, 
and  each  number  represents  an  electron;  thus  the  SGA  involves  permuting  electrons 
among  orbitals.  Flence  the  number  of  electrons  N  matches  the  number  of  boxes  in 


The  symmetrizer  and  antisymmetrizer 

Two  operators,  the  symmetrizer  and  the  antisymmetrizer,  relate  the  standard 
Young  Tableauxto  the  operations  they  label  [39],  The  symmetrizer  is  the  product  over 
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Figure  50.  Young  Tableaux  for  S5 


Figure  51.  Example  basis  vector  label  for  the  S5  irrep 


all  rows  of  the  sum  of  all  possible  permutations  in  a  single  row  of  a  tableau.  The 
antisymmetrizer  is  the  product  over  all  columns  of  the  sum  of  all  possible  even 
permutations  and  negative  odd  permutations  in  a  single  column  of  thetableau.  For 
example,  for  the  tableau 


the  first  row,  consisting  of  1  and  3,  includes  the  possible  permutations  (l)  and  (l  3) ; 
the  second  row  only  has  the  number  2,  and  so  the  only  permutation  is  the  identity  (2) . 
Thus  the  symmetrizer  is 

(2) [(1)+(13)]=(1)+(13)  (490) 

The  first  column,  consisting  of  1  and  2,  includes  the  possible  permutations  (l) ,  which  is 
even,  and  (l  2),  which  is  odd.  The  second  column  only  has  the  number  3,  so  the  only 
permutation  is  identity,  (3),  which  is  even.  Thus  the  antisymmetrizer  is 

(3) [(1)-(1  2)]  =  (1)-(1  2)  (491) 

The  operation  that  the  tableau  labels  is  the  product  of  the  symmetrizer  and  the 
antisymmetrizer.  Thus  the  above  tableau  labels  the  operation 

[W+f1  3)][(1)-(1  2)]  =  (1)-(1  2  3)  +  (l  3)-(l  2).  (492) 
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An  example:  The  S3  irreducible  representation 

Consider  as  a  thorough  example  the  group  S3 .  Using  the  above  prescription,  the 


Young  tableaux  are  associated  with  the  operations 


(1)  +  (1  2  3)  +  (l  3  2)  +  (2  3)  +  (l  3)  +  (l  2) 
(1)-(13  2)-(1  3)+(l  2) 

(l)-(l  2  3)+(l  3)  (I  2) 

(1)+(1  2  3)+(l  3  2)— (2  3) — (1  3) — (1  2) 


Note  that  all  these  operations  are  orthogonal  and,  to  within  a  constant,  idempotent. 
Using  the  formulation  in  equation  (488),  these  tableaux,  which  are  now  seen  to  be 
idempotent  projectors,  can  be  used  to  create  the  basis  vectors 


~  iK1  1  1  1  1  i) 

<->  1(1  -10  0-1  1) 


<->  1(1  0  -1  0  1  -1) 
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°  ±(‘  1  1  -1  -1  -1) 

where  normalization  has  been  enforced.  Note  that  there  are  only  four  basis  vectors 

here,  but  S3  has  6.  The  two  projectors  from  the  dog-leg  irrep  can  further  be  split  each 

into  a  pair  of  nilpotent  projectors.  Consider  the  operation 

1  I  2  |  (2  3) 

3 


which  yields  the  operation  proportional  to 


1  2  J^_3 

3  2 


_LJ]„(1  2  3)-(13  2)  +  (2  3) — (1  3) 


Likewise,  the  operation 


(493) 


yields  the  operation  proportional  to 


^-(1  2  3)  +  (l  3  2)  +  (2  3)-(l  2) 


(494) 


Both  of  these  operations  are  nilpotent  and  are  said  to  split  the  two  projectors  they  use 


because  they  have  the  relationship 
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These  projectors  can  be  used  to  construct  two  additional  basis  vectors,  just  like  the 
idempotent  projectors.  This  action  brings  the  total  number  of  basis  vectors  to  six,  as 
required.  These  six  vectors  will  block  diagonalize  all  permutations  of  S3  into  alxl 
symmetric  block,  a  lxl  antisymmetric  block,  and  two  2x2  blocks.  A  Hamiltonian 
which  commutes  will  all  those  permutations  will  also  be  thus  block  diagonalized  in  this 
basis.  The  interested  reader  may  wish  to  consult  Harter  [37]  for  a  parallel  approach  to 
finding  a  basis  forC3v. 

Spin  functions 

As  in  the  UGA,  in  the  SGA  each  eigenfunction  can  be  labeled  by  a  pair  of 
conjugate  tableaux,  one  spatial  and  the  other  spin.  Although  the  SGA  uses  the  Young 
tableaux  for  the  spatial  piece,  it  uses  the  Weyl  tableaux  for  the  spin  piece,  just  like  the 
UGA.  Figure  52  [41]  combines  the  information  from  Figure  8  and  Figure  49.  Recall  that 
in  the  UGA,  although  a  spin  state  could  be  identified  by  the  Young  frame,  the  specific 
path  taken  to  that  state  was  ambiguous.  In  the  SGA,  that  ambiguity  is  resolved. 
Surprisingly,  it  is  the  Young  tableau  of  the  spatial  state  rather  than  the  Weyl  tableau  of 
the  spin  state  that  specifies  this  path.  As  an  example,  consider  the  triplet  whose  irrep  is 
represented  by  the  frame 


and  thus  whose  conjugate  spatial  frame  is 
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Figure  52.  Genealogical  construction  of  spin  functions  labled  by  multiplicities  and 

Young  frames 


From  Figure  52  it  is  clear  that  this  spin  state  will  a  be  triply  degenerate  triplet.  While  the 
spin  state's  Weyl  tableau  determines  thez-projection  ms,  the  spatial  state's  Young 
tableau  determines  which  triplet  is  involved.  The  three  Young  tableaux. 
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each  specify  a  genealogical  path.  Take  the  first  tableau,  for  example.  1  is  in  the  first 
column,  so  we  add  it  to  the  total  spin.  2  appears  in  the  second  column,  so  we  subtract 
it.  3  and  4  are  both  found  in  the  first  column,  so  we  add  both  their  spins.  The  resultant 
path  of  add,  subtract,  add,  add  is  shown  in  Figure  53  [41].  Using  spin  coupling 
techniques,  we  recognize  that  this  coupling  path,  doublet,  singlet,  doublet,  triplet,  is  the 
triplet  of  spin  functions 

aj3aa)  -  \  J3aaa)^j 

<  l/2(\a/3a/3) -  \/3aa/3)  +  \a/3/3a) -  \/3a/3a)^  (495) 

yri(\am)~\PccPP)) 

The  second  tableau  then  represents  a  doublet,  triplet,  doublet,  triplet  path  outlined  in 
Figure  54  [41].  This  is  the  triplet 

|  aaj3a )  -  *Jy6  ( |  aj3aa )  + 1  f3aaa )) 

<  <Jy(jaaj3j3)-\j3j3aa)^  +  2J}h(\(xj3j3a)  +  \j3aj3a)-\aj3aj3)-\j3aaj3)^  (496) 

4)i(\am)+\PaW))-Jy,\WaP) 

Finally,  thethird  Young  tableau  represents  the  doublet,  triplet,  quartet,  triplet  coupling 
in  Figure  55  [41]  which  corresponds  to  the  triplet 

|  aaafi)  +  (|  aj3aa)  + 1  j3aao)  -  \  aafiot)} 

<  ‘Jy6{\[aaf30)  +  \af3af3)  +  \f3aaf3) -  \af3f3a) -  \f3a/3a) -\f3f3aa)^  (497) 

•M  ( I  »m) + 1  me)  - 1  me))  -  -M  (m« ) 

Although  all  three  triplets  are  clearly  different  sets  of  functions,  the  spin -orbit 
Flamiltonian  does  not  distinguish  between  them;  in  fact,  any  set  of  them  could  be  mixed 
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and  retain  a  definite  spin,  and  any  with  the  same  ms  mixed  would  also  retain  a  definite 
z-projection.  For  this  reason,  it  is  not  necessary  to  specify  exactly  which  genealogica  I 
path  was  taken  to  get  to  the  spin  state,  and  thus  why  the  UGA  does  not  suffer  despite 
this  disadvantage. 


electrons 


Figure  53.  Spin  state  represented  by  the  first  Young  Tableau 


electrons 
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Figure  54.  Spin  state  represented  by  the  second  Young  Tableau 


electrons 


Figure  55.  Spin  state  represented  by  the  third  Young  Tableau 
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